{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OLS Model practise - Data Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distances is minimum.  \n",
    "![img](https://pic2.zhimg.com/80/v2-d0ca37636d93b281c8eaee6ac4e4c4f1_720w.png)  \n",
    "\n",
    "The model is   \n",
    "![](https://pic3.zhimg.com/80/v2-31f810907af6fb3f87f3cc037cef174a_720w.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.          0.1010101   0.2020202   0.3030303   0.4040404   0.50505051\n",
      "  0.60606061  0.70707071  0.80808081  0.90909091  1.01010101  1.11111111\n",
      "  1.21212121  1.31313131  1.41414141  1.51515152  1.61616162  1.71717172\n",
      "  1.81818182  1.91919192  2.02020202  2.12121212  2.22222222  2.32323232\n",
      "  2.42424242  2.52525253  2.62626263  2.72727273  2.82828283  2.92929293\n",
      "  3.03030303  3.13131313  3.23232323  3.33333333  3.43434343  3.53535354\n",
      "  3.63636364  3.73737374  3.83838384  3.93939394  4.04040404  4.14141414\n",
      "  4.24242424  4.34343434  4.44444444  4.54545455  4.64646465  4.74747475\n",
      "  4.84848485  4.94949495  5.05050505  5.15151515  5.25252525  5.35353535\n",
      "  5.45454545  5.55555556  5.65656566  5.75757576  5.85858586  5.95959596\n",
      "  6.06060606  6.16161616  6.26262626  6.36363636  6.46464646  6.56565657\n",
      "  6.66666667  6.76767677  6.86868687  6.96969697  7.07070707  7.17171717\n",
      "  7.27272727  7.37373737  7.47474747  7.57575758  7.67676768  7.77777778\n",
      "  7.87878788  7.97979798  8.08080808  8.18181818  8.28282828  8.38383838\n",
      "  8.48484848  8.58585859  8.68686869  8.78787879  8.88888889  8.98989899\n",
      "  9.09090909  9.19191919  9.29292929  9.39393939  9.49494949  9.5959596\n",
      "  9.6969697   9.7979798   9.8989899  10.        ]\n"
     ]
    }
   ],
   "source": [
    "nsample = 100\n",
    "x = np.linspace(0,10,nsample) # Set array\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Constant item, set beta as beta_0, beta_1 to 1, 10\n",
    "X = sm.add_constant(x)\n",
    "beta = np.array([1,10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Error item - set length of normal distribution as k\n",
    "e = np.random.normal(size=nsample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "## y(t)\n",
    "y = np.dot(X, beta) + e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.00000</td>\n",
       "      <td>1.667889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.10101</td>\n",
       "      <td>0.967527</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.20202</td>\n",
       "      <td>2.018076</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.30303</td>\n",
       "      <td>2.360419</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.40404</td>\n",
       "      <td>3.848655</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         x         y\n",
       "0  0.00000  1.667889\n",
       "1  0.10101  0.967527\n",
       "2  0.20202  2.018076\n",
       "3  0.30303  2.360419\n",
       "4  0.40404  3.848655"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.DataFrame({\n",
    "    \"x\": x,\n",
    "    \"y\": y\n",
    "})\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.52193909 10.10222873]\n"
     ]
    }
   ],
   "source": [
    "model = sm.OLS(y,X)\n",
    "results = model.fit()\n",
    "print(results.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.999\n",
      "Model:                            OLS   Adj. R-squared:                  0.999\n",
      "Method:                 Least Squares   F-statistic:                 8.157e+04\n",
      "Date:                Sun, 26 Apr 2020   Prob (F-statistic):          6.08e-145\n",
      "Time:                        22:18:14   Log-Likelihood:                -143.97\n",
      "No. Observations:                 100   AIC:                             291.9\n",
      "Df Residuals:                      98   BIC:                             297.1\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.5219      0.205      2.549      0.012       0.116       0.928\n",
      "x1            10.1022      0.035    285.613      0.000      10.032      10.172\n",
      "==============================================================================\n",
      "Omnibus:                        2.218   Durbin-Watson:                   1.909\n",
      "Prob(Omnibus):                  0.330   Jarque-Bera (JB):                1.651\n",
      "Skew:                          -0.283   Prob(JB):                        0.438\n",
      "Kurtosis:                       3.274   Cond. No.                         11.7\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Find the fittedvalues get the y\n",
    "y_fitted = results.fittedvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.25, 2, -1, 20)"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAHWCAYAAADZ+5VfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXTU1d3H8fclBhgQExdUIGiVYtSyakQDiJGoEWUzKoobqC3uaPs0VlqLFhe0cWlR64bWpSpFiREUjTSIuKRqMEhQjIhaSVBBNCHIYIbkPn/cBAgkkGRm8puZfF7ncGbmzm9mvpzzPPjpvff3vcZai4iIiIh4p53XBYiIiIi0dQpkIiIiIh5TIBMRERHxmAKZiIiIiMcUyEREREQ8pkAmIiIi4rHdBjJjTE9jzBvGmE+MMR8bY66tHd/HGLPAGLOy9nHvRj4/ofaalcaYCaH+C4iIiIhEO7O7PmTGmG5AN2vth8aYLsASYCwwEfjBWnuHMeYGYG9r7R92+Ow+QCGQAtjazx5trf0x5H8TERERkSi12xkya+031toPa59XAiuAHsAY4Mnay57EhbQdZQALrLU/1IawBcCpoShcREREJFY0aw+ZMeYXwEDgPeAAa+03tW99CxzQwEd6AKu3e11aOyYiIiIitfZo6oXGmD2BOcB11toNxpit71lrrTEmqDOYjDGTgEkAnTt3Pvrwww8P5utEREREdvb997C6dq6ob1/Yo8lRqFFLliz53lrbNZjvaFIVxph4XBh7xlqbUzv8nTGmm7X2m9p9Zmsb+GgZkLbd6yRgUUO/Ya19BHgEICUlxRYWFjbpLyAiIiLSJLfeCn/+s3vesSM8+CCkpgb9tcaY/wX7HU25y9IAjwErrLX3bPfWXKDurskJwEsNfDwPOMUYs3ftXZin1I6JiIiItJ7p07eFMYBAABYt8qycHTVlD9kQ4EJguDFmae2f04A7gJONMSuBk2pfY4xJMcbMBLDW/gDcAnxQ+2da7ZiIiIhI69l7bzj1VPD5IC4O2reHtDSvq9pqt20vvKAlSxEREQna2rXw8cdw4onutbXw3/+6mbG0tJAsVwIYY5ZYa1OC+Y7gd7K1kkAgQGlpKZs3b/a6lLDp2LEjSUlJxMfHe12KiIhIdFu+HEaNgo0b4csvYc89wRgXwkIUxEIpagJZaWkpXbp04Re/+AXb3+EZK6y1rF+/ntLSUg455BCvyxEREYler70G48ZB584wf74LYxEuas6y3Lx5M/vuu29MhjEAYwz77rtvTM8AioiIhN1998Hpp0OvXvD++3DMMV5X1CRRE8iAmA1jdWL97yciIhJW1sKKFTByJLz1FvTs6XVFTRZVgSyS3Hzzzdx1112Nvp+bm8snn3zSihWJiIi0UeXl8Nlnbo/YjBmQkxMVy5Tbi9lAlltUxpA7FnLIDa8w5I6F5BaVte7vK5CJiIiE3xdfwODBbpkyEHCd9+PivK6q2WIykOUWlTElp5iycj8WKCv3MyWnOOhQdtttt3HYYYcxdOhQSkpKAHj00Uc55phj6N+/P2eeeSabNm3i3XffZe7cuWRlZTFgwABWrVrV4HUiIiIShLfegkGD4Ntv4dFHIYq7FMRkIMvOK8EfqK435g9Uk51X0uLvXLJkCbNmzWLp0qXMnz+fDz74AIDMzEw++OADPvroI4444ggee+wxBg8ezOjRo8nOzmbp0qX06tWrwetEJPZ5PVsvEhWWzYZ7+8DNie5x2ezdf+appyA9HfbdF957L6KavLZE1LS9aI415f5mjTfFW2+9xRlnnEGnTp0AGD16NADLly/nxhtvpLy8nI0bN5KRkdHg55t6nYjEjrrZ+rr/gVg3Ww8wdmAPL0sTiRzLZsO8yRCo/W90xWr3GqDfuIY/U1MD//wnDB0KL7wA++zTOrWGUUzOkHVP9DVrPBgTJ07k/vvvp7i4mJtuuqnRthVNvU5EYkc4ZutFYk7+tG1hrE7A78Z3tGkT/PADtGsHubmQlxcTYQxiNJBlZSTji6+/oc8XH0dWRnKLv3PYsGHk5ubi9/uprKxk3rx5AFRWVtKtWzcCgQDPPPPM1uu7dOlCZWXl1teNXScisSscs/UiMaeitGnja9bAsGFw5pmuvUVCQlTvGdtRTAaysQN7MD2zLz0SfRigR6KP6Zl9g1oiOOqoozjnnHPo378/I0aM4JjaRnO33HILxx57LEOGDOHwww/fev25555LdnY2AwcOZNWqVY1eJyKxqzVn60WiVkLS7seLitzm/U8/hd/9zrW3iDFRc7j4ihUrOOKIIzyqqPW0lb+nSFuw4x4ycLP1wf4PRJGYsuMeMoB4H4ya4faQvfQSnHee27w/bx707+9drY1oU4eLi4hEm7rQlZ1XwppyP90TfWRlJCuMiWyvbuN+/jS3TJmQBOlT3fjPP7sZsT59XDA78EBvaw0jBTIRkTAaO7CHApjI7vQbV/+Oyqoq1+S1QwdYsAC6dQNfbC/1x+QeMhEREYlS69fDKafANde414ceGvNhDBTIREREJBIUFLjlyf793fPjj/e6olalJUsRERHxVkEBnHii2zMG8PDDcP753tbUyjRDJiIiIt567bVtYSwuzi1btjEKZC108803c9dddzX6fm5uLp988kkrViQiIhJlamrc46mnug38cXHQvn3Un0vZErEbyFpyUGkIKZCJiIjsQmUljBkD//gHpKbCG2/ALbdAfr573cbEZiCrazJXsRqw2w4qDTKU3XbbbRx22GEMHTqUkhJ3Ft2jjz7KMcccQ//+/TnzzDPZtGkT7777LnPnziUrK4sBAwawatWqBq8TERFpk77+2h0M/uqr7lxKcCFsypQ2GcYgVgNZcw4qbaIlS5Ywa9Ysli5dyvz58/nggw8AyMzM5IMPPuCjjz7iiCOO4LHHHmPw4MGMHj2a7Oxsli5dSq9evRq8TkREpM157z13DNJXX8H8+XD55V5XFBFi8y7Lph5U2gxvvfUWZ5xxBp06dQJg9OjRACxfvpwbb7yR8vJyNm7cSEZGRoOfb+p1IiIiMeubb9zdlN26wcKFcOSRXlcUMWJzhqwpB5WGyMSJE7n//vspLi7mpptuYvPmzUFdJyIiErO6dYNHHnGzZApj9cRmIEuf6g4m3V68z4230LBhw8jNzcXv91NZWcm8efMAqKyspFu3bgQCAZ555pmt13fp0oXKysqtrxu7TkREJKZt3gwTJrhN+wAXXAD77edtTREoNgNZv3HulPiEnoBxj3WnxrfQUUcdxTnnnEP//v0ZMWIExxxzDAC33HILxx57LEOGDOHwww/fev25555LdnY2AwcOZNWqVY1eJyIiErO++w6GD4ennoKPPvK6mohmrLVe17CTlJQUW1hYWG9sxYoVHHHEER5V1Hrayt9TRERi3PLlMHIkrF3rAtlZZ3ldUdgYY5ZYa1OC+Y7Y3NQvIiIi3ikpgcGDYc89YfFiSAkqq7QJsblkKSIiIt457DC47jp4/32FsSZSIBMREZHgBQJw/fXwxRdgDEybBkmh724QqxTIREREJDjl5XD66ZCdDXPnel1NVNIeMhEREWm5Vavc5v1Vq+Dxx+Hii72uKCopkImIhFFuURnZeSWsKffTPdFHVkYyYwf28LoskaZbNtsdPVhR6hqsp0/d1kZq2TLX1sJaWLAATjjB21qjmJYsm6G0tJQxY8bQu3dvevXqxbXXXktVVRWLFi1i5MiRO13/8ssvM3DgQPr378+RRx7Jww8/7EHVIuKV3KIypuQUU1buxwJl5X6m5BSTW1TmdWkiTbNsNsybDBWrAese50124wC9esHJJ7vO+wpjQVEgayJrLZmZmYwdO5aVK1fy2WefsXHjRv70pz81eH0gEGDSpEnMmzePjz76iKKiItLS0lq3aBHxVHZeCf5Adb0xf6Ca7LwSjyoSaab8aRDw1x+r2gQ3XQuVldC5Mzz3HPzyl97UF0NiO5AVFMD06e4xSAsXLqRjx45cXLs2HhcXx7333svjjz/Opk2bdrq+srKSLVu2sO+++wLQoUMHkpOTg65DRKLHmnJ/s8ZFIk5Faf3XVRae90Put6BjAEMqeveQNTTbNG4cXHklbNoEQ4a4te2aGmjXDvr1g2uvhYkT4fvvd+4YvGjRLn/u448/5uijj643ttdee3HQQQfx+eef73T9Pvvsw+jRozn44INJT09n5MiRjB8/nnbtYjsDi8g23RN9lDUQvron+hq4WiQCJSS5ZcrVW6BkC3wagB8sjDkQLrvM6+piSuymg4oKF8bAPVZUtHoJM2fOJD8/n0GDBnHXXXdxySWXtHoNIuKdrIxkfPFx9cZ88XFkZWi2XKJE+lT4Jg6e3ATvVMF6Cxl7wrS/u15jEjLRO0O2qxmtTp3cVGp6OlRVQfv27nVqqnt/v/12OyO2oyOPPJIXXnih3tiGDRv4+uuv+eUvf8nrr7/e4Of69u1L3759ufDCCznkkEN44oknmvW7IhK96u6m1F2WErX6jYOOL0LNLPfaAIeN3HaXpYRM7M6QpaZCfj7ccot7rAtjLZSens6mTZt46qmnAKiurub//u//mDhxIp06ddrp+o0bN7Jou9C3dOlSDj744KBqEJHoM3ZgD965YThf3nE679wwXGFMooe1kJMD466Gjj6Ii3OP5072urKYFLuBDFwImzIl6DAGYIzhxRdf5Pnnn6d3794cdthhdOzYkdtvvx2A/Px8kpKStv4pKirir3/9K8nJyQwYMICbbrpJs2MiIhIdqqrgkkvgzDPh669DOsEhDYveJUsP9OzZk3nz5u00npaWht+/88bd448/vjXKEhERCZ3vv3dBbPFimDoVzjnH3RynIBZWuw1kxpjHgZHAWmttn9qxfwN1u1ITgXJr7YAGPvsVUAlUA1ustTryXUREJFJ9+qk7Bqm01O29Pu88rytqM5oyQ/YEcD/wVN2AtfacuufGmLuBXd3CeKK19vuWFigiIiKt5Pvv3XLlG29oRqyV7TaQWWsXG2N+0dB7xhgDjAOGh7YsERERaTVLlsDRR8PQobByJXTo4HVFbU6wm/qPB76z1q5s5H0LvG6MWWKMmRTkb2GtDfYrIlqs//1ERCTCVFfDdddBSorbsA8KYx4JdlP/eOC5Xbw/1FpbZozZH1hgjPnUWru4oQtrA9skgIMOOmin9zt27Mj69evZd999MTHYjM5ay/r16+nYsaPXpYiISFuwYQOMHw/z57uTbHQ4uKdaHMiMMXsAmcDRjV1jrS2rfVxrjHkRGAQ0GMistY8AjwCkpKTsNFWUlJREaWkp69ata2nJEa9jx44kJSV5XYaIiMS6r76CUaNgxQp48EG4/HKvK2rzgpkhOwn41Fpb2tCbxpjOQDtrbWXt81OAaS39sfj4eA455JCWflxERETqvPkmrF4Nr74KJ5/sdTVCE/aQGWOeAwqAZGNMqTHm0tq3zmWH5UpjTHdjzPzalwcAbxtjPgLeB16x1r4WutJFRCJfblEZQ+5YyCE3vMKQOxaSW1TmdUkSS5bNhnv7wM2J7nHZ7F1fv2aNe5wwAT77TGEsgphI3EiekpJiCwsLvS5DRCQouUVlTMkpxh+o3jrmi49jemZfHaEkwVs2G+ZNhsB2jcnjfTBqxs5nTVoL06bBX/8K770Hffq0bq0xzhizJNheq7F9dJKIiIey80rqhTEAf6Ca7LwSjyqSmJI/rX4YA/c6f4fdQZs3w/nnw803w7hx0Lt3q5UoTaejk0REwmRN+c5Hqu1qXKRZKhrcwl1//LvvYOxYNyt2552QlQUx2KkgFmiGTEQkTLon+po1LtIsCY3clb/9+D/+AcuWwZw5cP31CmMRTIFMRCRMsjKS8cXH1RvzxceRlZHcyCdEmiF9qtsztr14nxvftMm9/vOfobAQzjij9euTZlEgExEJk7EDezA9sy89En0YoEeiTxv6JXT6jXMb+BN6AsY9jvw7LFwDv/oVfPst7LEHHHGE15VKE2gPmYhIGI0d2EMBTMKn37htd1QGAnDNNfDww5CZCV26eFubNIsCmYiISDQrKHANXl991S1PTpkCt94K7bQIFk0UyERERKJVQQGkp7vWFtbCjTfCLbd4XZW0gOKziIhItHrjDaiqcmGsXTvo1MnriqSFFMhERESi0T//CTk50L49xMVBhw6QluZ1VdJCWrIUERGJJjU1bp/YX/8KJ50EL73k9o6lpUFqqtfVSQspkImIiESLn36CCy6A3Fy47DK47z6Ij9ch4TFAS5YiIiLRYsIEmDsX/vY3ePBBF8YkJmiGTEREJFrccgtccgmcdprXlUiIaYZMREQkkuXkuIav1rqu+wpjMUmBTEREJBJZC9Onw5lnuk37dedTSkxSIBMREYk0P/8MEyfCH/8I48e7fmOdO3tdlYSRApmIiEikycyEp56Cv/wFnnkGOnb0uiIJM23qFxERiTRXXgkXXgjnnut1JdJKFMhEREQiwYIF8OWXMGkSnH6619VIK9OSpYiIiNcefBBGjICHHoJAwOtqxAMKZCIiIl7ZsgWuvdYtUZ56Krz5ppq9tlFashQREfFCTQ2MHQuvvAK//S1kZ7tDwqVNUiATERHxQrt27kDwUaPcuZTSpimQiYiItKaCAvD7Yfhw+P3vva5GIoT2kImIiLSWZ5+FE0+EG25wnfhFammGTEQkjHKLysjOK2FNuZ/uiT6yMpIZO7CH12VJa1g2G/KnQUUp7NUDPjsCHp4DJ5wAc+aAMV5XKBFEgUxEJExyi8qYklOMP1ANQFm5nyk5xQAKZbFu2WyYNxm+qIQvt8BXJfDlJzAmDWbnQfv2XlcoEUZLliIiYZKdV7I1jNXxB6rJzivxqCJpNfnTXBh7ahMsqoKvqiElHoatVRiTBimQiYiEyZpyf7PGJYZUlMLyAFQDdVvF9moHG8q8rEoimAKZiEiYdE/0NWtcYkhpAiyp7bhvgDjgF3GQkORlVRLBFMhERMIkKyMZX3z9Rp+++DiyMpI9qkjCzlq45x54fDUcsAeM88GJHeCiTnBoF0if6nWFEqG0qV9EJEzqNu7rLss2IhCAq66CRx+FzEzIGgMFf3XLlwlJLoz1G+d1lRKhjI3APigpKSm2sLDQ6zJERESabtMm13n/pJPg1ltdJ35pE4wxS6y1KcF8h2bIREREgvH557D//rDXXrB4MXTs6HVFEoUU30VERFpq0SIYNAiuvNK9VhiTFlIgExERaYnHH4eTT4YDD4S//MXraiTKKZCJiIg0R00NXH89XHqpO5fy3XehVy+vq5Iop0AmIiLSHOvWwdNPu2XK+fMhMdHriiQGaFO/iIhIU3z3HXTtCgccAEuXukeRENEMmYiIyO4UFsLAgXDTTe61wpiEmAKZiIjIrsyZA8OGQYcOcM45XlcjMUqBTEREpCHWwu23w1lnudmx996DPn28rkpi1G4DmTHmcWPMWmPM8u3GbjbGlBljltb+Oa2Rz55qjCkxxnxujLkhlIWLiIiE1cqVrp3F+edDfr5r/ioSJk2ZIXsCOLWB8XuttQNq/8zf8U1jTBzwADACOBIYb4w5MphiRUREws7vd4+HHeb2jj39tBq+StjtNpBZaxcDP7TguwcBn1trv7DWVgGzgDEt+B4REZHW8cknblnyX/9yr/v2BWO8rUnahGD2kF1tjFlWu6S5dwPv9wBWb/e6tHasQcaYScaYQmNM4bp164IoS0REpAXy8iA1FX76yc2OibSilgayB4FewADgG+DuYAux1j5irU2x1qZ07do12K8TERFpugcegNNPh1/8At5/351PKdKKWtQY1lr7Xd1zY8yjwMsNXFYG9NzudVLtmIiIiPcKCtzh4PvtB1dfDSNHwrPPQpcuXlcmbVCLApkxppu19pval2cAyxu47AOgtzHmEFwQOxc4r0VVioiIhFJBAaSnQ1UVtG8P2dnw299CXJzXlUkbtdtAZox5DkgD9jPGlAI3AWnGmAGABb4CLqu9tjsw01p7mrV2izHmaiAPiAMet9Z+HJa/hYiISHPk5Gy7m7KqCgIBhTHx1G4DmbV2fAPDjzVy7RrgtO1ezwd2aokhIiLimXfegUcfdc/btXMzZGlpnpYkosPFRUSk7fjXv+DSS+Ggg+Dhh+GLL1wYS031ujJp4xTIRESkbZg/Hy680AWwF16Afff1uiKRrXSWpYiItA0ZGTBjhus3pjAmEUaBTEREYte330JmJqxZ4zbtX3ON2zMmEmG0ZCkiEka5RWVk55WwptxP90QfWRnJjB3Y6KElEkoffQSjRsEPP8Cnn0L37tveWzYb8qdBRSkkJEH6VOg3zrtapc1TIBMRCZPcojKm5BTjD1QDUFbuZ0pOMYBCWbjNmwfjx8Pee8Pbb8OAAdveWzYb5k2GQG3bi4rV7jUolIlntGQpIhIm2XklW8NYHX+gmuy8Eo8qaiPmzIExY+DII90xSNuHMXAzY3VhrE7A78ZFPKJAJiISJmvK/c0alxBJT4ff/c4di9St287vV5Q2/LnGxkVagQKZiEiYdE/0NWtcgvDDDy6Ebd4MiYlw113QqVPD1yYkNW9cpBUokImIhElWRjK++PrH8fji48jKSPaoohj12Wdw3HHwwANuiXJ30qdC/A6hON7nxkU8ok39IiJhUrdxX3dZhtEbb8CZZ7qWFgsXwpAhu/9M3cZ93WUpEcRYa72uYScpKSm2sLDQ6zJERCSSzZ4N558Phx3m7qo89FCvK5I2yhizxFqbEsx3aMlSRESiU//+bnbs3XcVxiTqKZCJiEj02LgR7r8frIXkZJg1CxISvK5KJGgKZCIiEh1Wr4ahQ+Haa+HDD72uRiSkFMhERCTyffABDBoEX34Jr7wCRx/tdUUiIaVAJiIikS0nB4YNg44d3X6xU0/1uiKRkFMgExGRyNapk5sde/99+NWvvK5GJCwUyEREJPK8+SZcfDEUFLgZsUWLoGtXr6sSCRsFMhERiSyvvALDh8MTT7jHggIwxuuqRMJKnfpFRMIot6gstjr1L5sd3g73H38MF1wANTXudSDgZsdSU0P3GyIRSIFMRCRMcovKmJJTjD9QDUBZuZ8pOcUA0RnKls2GeZMh4HevK1a71xCaUPbaazBuHMTHQ4cOsGULtG8PaWnBf7dIhNOSpYhImGTnlWwNY3X8gWqy80o8qihI+dO2hbE6Ab8bD4VPPnEd95cudWdU3nIL5OdrdkzaBM2QiYiEyZpyf7PGI15FafPGm2LLFigpcXdP/va3cOWVrr1Fz54KYtKmaIZMRCRMuif6mjUe8RKSmje+OxUVcPrpMHgwrF3rNu537OiWRu/tAzcnusdls1tes0iUUCATEQmTEw9vuE1DY+MRL30qxO8QJuN9bry5vvjCBbGFC+Gee2D//d143T61itWA3bZPTaFMYpwCmYhImLzx6bpmjUe8fuNg1AxI6AkY9zhqRvM39L/9Nhx7LHzzDbz+Olx66bb3wr1PTSRCaQ+ZiEiYxNweMnDhK9g7Kp94Avbe2/Ub6927/nvh2KcmEgU0QyYiEiYxt4csGDU1bp8YwAMPwHvv7RzGIPT71ESihAKZiEiYZGUk44uPqzfmi48jKyPZo4o8smkTnHMODB0KGze6HmN7793wtaHcpyYSRbRkKSISJnXNX2OqU39zffMNjB4NS5bAXXdB5867vr5uOTScpwGIRCBjrfW6hp2kpKTYwsJCr8sQEZFgLF0Ko0bBjz/Cc8+55yIxyBizxFqbEsx3aIZMRETC4/e/d73F3nkH+vf3uhqRiKZAJiIioWMt/Pyza/D6zDNuM3+3bl5XJRLxtKlfRERCo6oKfv1rGDPGHYl0wAEKYyJNpEAmIiLBW78eTjkFHn8cBg2CdvrPi0hzaMlSRESCU1ICI0fC11/D00/DBRd4XZFI1FEgExGRlqupgbPOcgeFL1wIQ4Z4XZFIVFIgExGRlrHWLU0+/TQkJMAhh3hdkUjU0iK/iIg0z9tvu677553nXg8YoDAmEiTNkImISNP95z+QkeGWKvfYA665BgYP9roqkainGTIREWmar7+GCy90YQzckuWbb3pbk0iMUCATEZHd+/lnOOEE2LDBHQ4eFwft20NamteVicSE3S5ZGmMeB0YCa621fWrHsoFRQBWwCrjYWlvewGe/AiqBamBLsOc8iYiIRzp0gL/9DXr3dndULlrkwlhqqteVicSE3R4ubowZBmwEntoukJ0CLLTWbjHG3Algrf1DA5/9Ckix1n7fnKJ0uLiISASwFm69FQ46CCZM8LoakYgVisPFd7tkaa1dDPyww9jr1tottS//CyQFU4SIiESYzZvdfrGpU91dlSISVqHYQ3YJ8Goj71ngdWPMEmPMpF19iTFmkjGm0BhTuG7duhCUJSIiLbJ2LQwf7g4Hv/12eOQRrysSiXlBtb0wxvwJ2AI808glQ621ZcaY/YEFxphPa2fcdmKtfQR4BNySZTB1iYhEityiMrLzSlhT7qd7oo+sjGTGDuwRnh9bNhvyp0FFKSQkQfpU6Deued+xYQMceyx89x288AKceWZ4ahWRelocyIwxE3Gb/dNtIxvRrLVltY9rjTEvAoOABgOZiEisyS0qY0pOMf5ANQBl5X6m5BQDhD6ULZsN8yZDwO9eV6x2r6F5oWyvveDKK+HEEyFF92GJtJYWLVkaY04FrgdGW2s3NXJNZ2NMl7rnwCnA8pYWKiISbbLzSraGsTr+QDXZeSWh/7H8advCWJ2A343vjrVw//1QUOBeZ2UpjIm0st0GMmPMc0ABkGyMKTXGXArcD3TBLUMuNcY8VHttd2PM/NqPHgC8bYz5CHgfeMVa+1pY/hYiIhFoTbm/WeNBqSht3nidQACuusp13H/88dDXJSJNstslS2vt+AaGH2vk2jXAabXPvwD6B1WdiEgU657oo6yB8NU90Rf6H0tIcsuUDY03prwcxo2DBQvg+uth+vTQ1yUiTaJO/SIiYZKVkYwvPq7emC8+jqyM5ND/WPpUiN8h6MX73HhD1q51Z1C+8QbMnAl33gnt9J8EEa/ocHERkTCp27jfKndZ1m3cb+pdlvvuC4MGwT/+oeOPRCLAbjv1e0Gd+kVEwuS55+D44yFJ/bxFQqVVOnw9Sz4AACAASURBVPWLiEgMqKmBKVPgvPPgr3/1uhoR2YGWLEVEYt1PP8FFF0FODkyaBHff7XVFIrIDBTIRkVj27bcwciR8+CHccw9cdx0Y43VVIrIDBTIRkVjm88Eee8BLL8GoUV5XIyKNUCATEYlF//kPDBkCCQmuA79mxUQimjb1i4jEknffhREj4OST4Y473JjCmEjE0wyZiEisWLwYhg+H6mqIi3MHhItIVNAMmYhILFi/HiZOdGGsTt1h4SIS8RTIRERiwQ8/QGUltG/vZsfat1cHfpEooiVLEZFoVlwMffpA797w9dewdCksWuTCWGqq19WJSBNphkxEJFo9/DAMHOgewbW4SE11HfkVxkSiigKZiEi0qa6G3/4WLr8cTjnFHYckIlFNgUxEJJpUVsKYMfC3v8G118LcubDXXl5XJSJB0h4yEZFoUlQE+fnwj3/AFVd4XY2IhIgCmYhINFi7FvbfH4YNgy++gG7dvK5IREJIS5YiIpFu1iw45BB49VX3WmFMJOYokImIRCpr4S9/gfHj4eij4ZhjvK5IRMJES5YiIpFo82a45BJ47jmYMMG1tujQweuqRCRMNEMmIhKJcnJcGJs+Hf75T4UxkRinGTIRkUjy888ufI0fD4cfDkcd5XVFItIKNEMmIhIpXnkFevWC5cvBGIUxkTZEgUxExGvWwt//DqNHu9YWiYleVyQirUyBTETES4EAXHklXHedC2RvvQVJSV5XJSKtTIFMRMRL990HDz0Ef/gDzJkDnTt7XZGIeECb+kVEvGCt2yd29dXQuzeMGuV1RSLiIc2QiYi0tsWLYfBgWL8e2rdXGBMRBTIRkVZTUABnnQXDh8OPP8KGDV5XJCIRQkuWIiJhlFtURnZeCWd/8k+ufeZfmBrAANdluvMpRUTQDJmISNjkFpUxJaeYozcs4Jq3ZrkwVifnPlg227PaRCSyKJCJiIRJdl4J/kA11+8xm7jB8e5fXAPEAT1rIH+axxWKSKTQkqWISJjs/Wkxvyt8ie5j10HveJho4Ktq+EUc9NwDKkq9LlFEIoQCmYhIOLz4Is8/ewPrfV34bkMi3fapcCGs53b/7CaoAayIOFqyFBEJJWvhzjshMxP/4Udy7iV/Z/qe57PJtq9/XbwP0qd6U6OIRBwFMhGRUJoyBW64AcaPZ5/33+H3Fw1jyV4nMyXwa76lKxYDCT1h1AzoN87rakUkQhhrrdc17CQlJcUWFhZ6XYaISPMtXQqvvAJ//KPrxC8iMc8Ys8RamxLMd2iGTEQkWCtWwPTp7vmAAfCnPymMiUizKJCJiARjwQJITYW//Q2++87rakQkSimQiYi01EMPwYgR0LMnvP8+HHCA1xWJSJRSIBMRaYkbboArroCMDHjnHTj4YK8rEpEo1qRAZox53Biz1hizfLuxfYwxC4wxK2sf927ksxNqr1lpjJkQqsJFRDyVkgLXXQdz58Jee3ldjYhEuabOkD0BnLrD2A1AvrW2N5Bf+7oeY8w+wE3AscAg4KbGgpuISMT76iuYM8c9P+ssuPdeiIvztCQRiQ1NCmTW2sXADzsMjwGerH3+JDC2gY9mAAustT9Ya38EFrBzsBMRiXwFBTBokFum3LjR62pEJMYEc3TSAdbab2qffws0tJu1B7B6u9eltWMiIk2zbLY7hLui1B01lD619RuqPvssXHIJJCXByy/Dnns2+aO5RWVk55WwptxP90QfWRnJjB2ofwZFpL6QbOq3rrtsUB1mjTGTjDGFxpjCdevWhaIsEYl2y2bDvMlQsRqw7nHeZDfeWm6+Gc4/H449Ft57Dw4/vMkfzS0qY0pOMWXlfixQVu5nSk4xuUVlYStXRKJTMIHsO2NMN4Dax7UNXFMG9NzudVLt2E6stY9Ya1OstSldu3YNoiwRiRn50yDgrz8W8Lvx1lJTAxMnun5j++7brI9m55XgD1TXG/MHqsnOKwlhgSISC4IJZHOBursmJwAvNXBNHnCKMWbv2s38p9SOiYjsXkVp88ZD5dtvYckS9/wvf4HHH4f27Xf9mQasKfc3a1xE2q6mtr14DigAko0xpcaYS4E7gJONMSuBk2pfY4xJMcbMBLDW/gDcAnxQ+2da7ZiIyO4lJDVvPBSWLXPLk2edBYGAOwKphccgdU/0NWtcRNqupt5lOd5a281aG2+tTbLWPmatXW+tTbfW9rbWnlQXtKy1hdbaX2/32cettb+s/fPPcP1FRCQGpU+F+B3CS7zPjYdaQQFMmADHHQdbtrj2FvHxQX1lVkYyvvj6bTF88XFkZSQH9b0iEnuCuctSRCS86u6mDPddlu++C2lp22bEHnwQjjoq6K+tu5tSd1mKyO4okIlIZOs3LvxtLhYtcmEMoF07+PhjGD06JF89dmAPBTAR2S0FMhFpu378ETZtghNPBJ8Pqqrc5v20NK8rE5E2RoFMRNqmzz+H00+HhATXXyw/382UpaVBaqrX1YlIG6NAJiJtz5tvQmam2y82c6Z7TE1VEBMRz4SkU7+ISNT45z/h5JNh//3dzNjxx3tdkYiIApmItCFVVXDPPXDCCa7NRa9eXlckIgJoyVJE2oKffnLLkp06wX/+A/vsE3SPMRGRUNIMmYjEttJStyx58cXu9QEHKIyJSMRRIBOR2FVYCIMGuTsqJ0zY/fUiIh5RIBOR2DRnDgwb5vqKvfMOnHaa1xWJiDRKgUxEYk9lJVxxBfTv7+6k7NvX64pERHZJm/pFJHZUVbn9YV26wMKF8MtfQseOXlclIrJbmiETkdjw/feQng7TprnXffoojIlI1FAgE5Hot2IFHHus28R/+OFeVyMi0mxashSR6LZgAZx9tpsNW7TIBTMRkSijGTIRiV7ffgujR8PBB8P77yuMiUjU0gyZiEQfa13n/QMPhJwcGDrUbeQXEYlSmiETkeiyYQOMGgUvvuhejxihMCYiUU+BTESix5dfwuDBkJcHP/7odTUiIiGjJUsRiXwFBfDkk/Dvf7vXr73mWlyIiMQIBTIRiWwFBTB8OGze7PaNPfuswpiIxBwtWYpIZFu0CAIB97xdO7dsKSISYxTIRCQy+f1wySWQlOQOCI+Lc49paV5XJiISclqyFJHI8+23MGYMfPABHHcc5Oe7mbK0NEhN9bo6EZGQUyATkcjy0UeurcX69a7H2NixblxBTERimAKZiESOoiIYNgwSEuDtt2HgQK8rEhFpFdpDJiKR41e/gosvdscgKYyJSBuiQCYi3goE4MYbYd06t2l/xgzo3t3rqkREWpUCmYh454cf4NRT4bbb4KWXvK5GRMQz2kMmEquWzYb8aVBRCglJkD4V+o3zuqptVq6EkSNdX7Enn4SLLvK6orDILSojO6+ENeV+uif6yMpIZuzAHl6XJSIRRoFMJBYtmw3zJkPA715XrHavITJC2fvvu5mxdu1cS4vjj/e6orDILSpjSk4x/kA1AGXlfqbkFAMolIlIPVqyFIlF+dO2hbE6Ab8bjwSHHuoOCX/vvZgNYwDZeSVbw1gdf6Ca7LwSjyoSkUilQCYSiypKmzfeGqqr4aGH3Cb+/faDl1+GXr28q2d7y2bDvX3g5kT3uGx2SL52Tbm/WeMi0nYpkInEooSk5o2H28aNkJkJV1wBL7zgTQ2NqVverVgN2G3LuyEIZQm++GaNi0jbpUAmEovSp0K8r/5YvM+Nt7bVq2HoUDcjdt99MH5869ewK2Fc3jWmeeMi0nZpU79ILKrbuO/1XZZLlrg7KX/6CV55xW3kjzRhXN4t3xRo1riItF0KZCKxqt847++ojI+Hrl1hwQLo08fbWhqTkFS7XNnAeJC6J/ooa2C/WPdEXwNXi0hbpiVLEQktayEvzz3v1w+WLo3cMAZhXd7NykjGFx9Xb8wXH0dWRnLQ3y0isUWBTERC5+efYcIEtzRZF8raRfg/M/3GwagZkNATMO5x1IyQzC6OHdiD6Zl96ZHowwA9En1Mz+yrHmQishNjrfW6hp2kpKTYwsJCr8sQkeZYtw7OOAPeeQduvRX++EftXheRNsEYs8RamxLMd2gPmYgEp6AAZs+Gf/8bfvzRPT/7bK+rEhGJKgpkItJyBQWQnu6WKmtqYOZMhTERkRZo8eYOY0yyMWbpdn82GGOu2+GaNGNMxXbXeNAESUTCZvZsqKpyYSwuDtau9boiEZGo1OIZMmttCTAAwBgTB5QBLzZw6VvW2pEt/R0RiUBbtsBvfwsPPgh71P4z0r49pKV5WpaISLQK1ZJlOrDKWvu/EH2fiESqigo45xx3F+VFI2Gv5fDJN9DnQOi8Gkj1usKIkltURnZeCWvK/XRP9JGVkay7LEVkJ6EKZOcCzzXyXqox5iNgDfB7a+3HIfpNEWltX37pOu9/9hlMnQTtX3LHDB3fAVjvzoAE7xvSRojcojKm5BTjD1QDUFbuZ0pOMYBCmYjUE3SDIGNMe2A08HwDb38IHGyt7Q/cB+Tu4nsmGWMKjTGF69atC7YsEQmH55+HNWvc7FjiO2E7AzJWZOeVbA1jdfyBarLzSjyqSEQiVSg6No4APrTWfrfjG9baDdbajbXP5wPxxpj9GvoSa+0j1toUa21K165dQ1CWiITM99+7x6wsKC6G4cPDegZkrFjTwLFJuxoXkbYrFIFsPI0sVxpjDjTGdYY0xgyq/b31IfhNEWkNNTVw441wxBHwv/+5Rq9JtWc8NnbWYwjOgIwVjZ1ZqbMsRWRHQQUyY0xn4GQgZ7uxy40xl9e+PAtYXruHbAZwro3EowFEZGebNrnN+7fdBmPGQLdu9d8P4xmQsUJnWYpIUwW1qd9a+xOw7w5jD233/H7g/mB+Q0Q88M03MHo0LFkCd90Fv/vdzscg1W3cz5/mlikTklwY04b+reo27usuSxHZHZ1lKSI7u+oqePJJePZZF8xERKRRoTjLMhR7yEQkVlRVucfsbPjvfxXGRERaiQKZiIC1cPfdMGgQbNgAnTpBnz5eVyUi0mYokIm0dVVVMGkS/P73cNhh245CEhGRVqNAJtKW/fADnHoqzJwJf/4zzJrlZsdERKRV6X8Ki7RlkybBO+/Av/4F55/vdTUiIm2WZshE2qK6u6vvvRfeeENhTETEYwpkIm3NzJmu4WtNDfTsCYMHe12RiEibp0Am0lZUV7uN+7/5jbuT0q/zFEVEIoUCmUhbsHEjnHGGa21xzTXw8svQubPXVYmISC1t6heJZQUFsGgRvPgifPgh3H+/68IvIiIRRYFMJFYVFEB6uusztscervu+wpiISERSIBOJVTNmwObN2+6o3LzZ23pERKRRCmQiscZauO021+S1XTv3p317SEvzujIREWmEAplILNm8GX79a3jmGbjgAvf83XddGEtN9bo6ERFphAKZSKyoqYERI9wm/ltvhT/+EYyBE07wujIREdkNBTKRWNGuHUycCFdeCWef7XU1IiLSDApkItHutdfgp5/gzDNhwgSvqxERkRZQY1iRaGUt3HcfnH66a/haU+N1RSIi0kIKZCLRaMsWuPpqmDwZRo6E1193S5YiIhKVtGQpEm2qqmDUKBfCsrJg+nSIi/O6KhERCYICmUi0ad8e+vWDcePg0ku9rkZEREJAgUwkWrz9Nuy5JwwY4I5BEhGRmKFNJyLR4Kmn3LmUv/+915WIiEgYKJCJRLKaGvjTn1w7i6FD4fnnva5IRETCQEuWIpHK74cLL4Q5c+A3v4EHHoD4eK+rEhGRMNAMmUik2mMP2LgR7rkHHn5YYUxEJIZphkwk0ixdCt27w/77w/z56i8mItIG6F96kUiSmwtDhrimr6AwJiLSRuhfe5FIYK1rZZGZCX36wIwZXlckIiKtSIFMxGtVVfDrX8P118PZZ8OiRXDggV5XJSIirUiBTMRrGzfC4sUwdSo89xz4fF5XJCIirUyb+kW88sIL8OmnruFrUZHrwi8iIm2SZshEvHDffW55cupUF8iKi72uSEREPKRAJtLaHnkErr3WPbfW7SFbtMjTkkRExFsKZCKtpboafvc7uOwyGDTI7RWLi4P27SEtzevqRETEQ9pDJtJavvoKZs6EyZPh7rvhgw/czFhaGqSmelyciIh4SYFMJNzWr4d99oFevWD5cjjoIDeemqogJiIigJYsRcLrvffgV79yB4PDtjAmIiKyHQUykXD597/hhBOgc2cYPtzrakREJIIpkImEmrUwbRqcey4cc4ybJTvySK+rEhGRCKZAJhJqhYVw881w0UXwn//Afvt5XZGIiEQ4beoXCZVAAOLj3axYQYFrbWGM11WJiEgUCHqGzBjzlTGm2Biz1BhT2MD7xhgzwxjzuTFmmTHmqGB/UyTiLF8ORxwBCxa418ceqzAmIiJNFqoZshOttd838t4IoHftn2OBB2sfRWLD/Pluv1iXLq69hYiISDO1xh6yMcBT1vkvkGiM6dYKvysSXtbCjBkwahT88pfw/vtw9NFeVyUiIlEoFIHMAq8bY5YYYyY18H4PYPV2r0trx0Si22uvuTMpR4+Gt96CHvo/axERaZlQLFkOtdaWGWP2BxYYYz611i5u7pfUhrlJAAepeaZEMmvd/rBTT4VZs+Dss6FdM/+3zbLZkD8NKkohIQnSp0K/ceGpVzyVW1RGdl4Ja8r9dE/0kZWRzNiBCu8iUl/QM2TW2rLax7XAi8CgHS4pA3pu9zqpdmzH73nEWptirU3p2rVrsGWJhMeqVXD88bBypQtl55zTsjA2bzJUrAase5w32Y1LTMktKmNKTjFl5X4sUFbuZ0pOMblFO/0TKCJtXFCBzBjT2RjTpe45cAqwfIfL5gIX1d5teRxQYa39JpjfFfHE4sXu7skVK2Dt2pZ/T/40CPjrjwX8blxiSnZeCf5Adb0xf6Ca7LwSjyoSkUgV7JLlAcCLxt3evwfwrLX2NWPM5QDW2oeA+cBpwOfAJuDiIH9TpPU98QRMmgSHHgovv+w28bdURWnzxiVqrSn3N2tcRNquoAKZtfYLoH8D4w9t99wCVwXzOyKemjULLr4Y0tPh+edh772D+76EpNrlygbGJaZ0T/RR1kD46p7o86AaEYlkOjpJZHfGjIE77oBXXw0+jIHbwB+/w3+Q431uXGJKVkYyvvi4emO++DiyMpI9qkhEIpUCmUhD5s6FAQNc532fD/7wB3csUij0GwejZkBCT8C4x1EzdJdlDBo7sAfTM/vSI9GHAXok+pie2Vd3WYrIToxbUYwsKSkptrBwp1OYRFrHE0/AJZe49hYdOsAbb0BqqtdViYhIhDLGLLHWpgTzHZohE9neiy+6zft1/0NlyxZYtMjTkkREJPYpkInUeeopyMyE3r2hY0eIi4P27SEtzevKREQkxoXqcHGR6JeeDlddBdnZsHSpmxlLS9NypYiIhJ32kEnbtn69OyB86lQ3IyYiItJM2kMmEoxPP3Wd9++8Ez76yOtqRESkDVMgk7bpP/+B446Dykq3NHnUUV5XJCIibZgCmbQ9Tz0Fp54KPXvC+++7YCYiIuIhBTJpe3r3dt3333kHDj7Y62pEREQUyKSN2LAB/vUv9zw1FebMgb328rYmERGRWgpkEvv+9z8YMsQdEL5ypdfViIiI7ESBTGLbf/8LgwbB6tXucPDevb2uSEREZCcKZBK7Zs1yjV27dHHB7KSTvK5IRESkQerUL7Hr559dn7E5c2C//Zr+uWWzIX8aVJRCQhKkT4V+48JXp8S03KIysvNKWFPup3uij6yMZMYO7OF1WSISYTRDJrHF74e333bPJ0yAhQubH8bmTYaK1YB1j/Mmu3GRZsotKmNKTjFl5X4sUFbuZ0pOMblFZV6XJiIRRoFMYsd338Hw4XDKKe45NP84pPxpEPDXHwv43bhIM2XnleAPVNcb8weqyc4r8agiEYlUWrKU2FBcDCNHwrp1rr3FAQe07HsqSps3LrILa8r9zRoXkbZLM2QS/V55BQYPhi1b4K23IDOz5d+VkNS8cZFd6J7oa9a4iLRdCmQS/d54Aw47zB2DdPTRwX1X+lSI3+E/lvE+Ny7STFkZyfji6y+b++LjyMpI9qgiEYlUWrKU6BQIuN5ihx4Kd97p7qjs1Cn47627m1J3WUoI1N1NqbssRWR3jLXW6xp2kpKSYgsLC70uQyJVXh5cdRVUVMAXX7g+YyIiIh4xxiyx1qYE8x1aspTo8vzzMGIErFrlzqdcvtzrikRERIKmQCbR4803XW+xulnd6mpYtMjTkkREREJBgUyix913Q9eu0LGj6y/Wvr07GklERCTKaVO/RLaaGrc0mZgITz/tXn/6qZsZS0uD1FSvKxQREQmaAplErp9+ggsugG++gcWLISHBjaemKoiJiEhM0ZKlRKayMjj+eJg7F847D+Ljva5IREQkbDRDJpFnyRIYPRoqK2HePDjtNK8rEhERCSsFMoksNTXw61+7GbF33oG+fb2uSEREJOwUyCQyWOvOooyPhxdegD33bPkB4SIiIlFGgUy89/PPcNllLpQ98QT06uV1RSIiIq1Km/rFW99/DyefDE8+qSAmIiJtlmbIxDsrVsDIkbBmDcyaBeec43VFIiIinlAgE29UVbkzKTdvdk1ejz3W64pEREQ8o0Amrc9ad+zRU0/BwQe7PyIiIm2Y9pBJ69myBa69Fu69170eNkxhTEREBAUyaS0bNrhmrzNmuC78IiIispWWLCX8vvrKbd7/9FN46CHX4kJERES2UiCT8KqsdAeB+/3w2mtw0kleVyQiIhJxFMgkvLp0gdtvd6Hs8MO9rkZERCQiaQ+ZhF5NDdx0k5sRA7j4YoUxERGRXWhxIDPG9DTGvGGM+cQY87Ex5toGrkkzxlQYY5bW/pkaXLkS8RYtggEDYNo0ePVVr6sRERGJCsEsWW4B/s9a+6ExpguwxBizwFr7yQ7XvWWtHRnE70i0ePllGDPGzZDFx6vzvoiISBO1eIbMWvuNtfbD2ueVwAqgR6gKkyjz3XdwwQUujIF7fPNNb2sSERGJEiHZQ2aM+QUwEHivgbdTjTEfGWNeNcb8ahffMckYU2iMKVy3bl0oypLWtP/+cPrp0KEDxMW5TvxpaaH/nWWz4d4+cHOie1w2O/S/IRJCuUVlDLljIYfc8ApD7lhIbpH68InIzoy1NrgvMGZP4E3gNmttzg7v7QXUWGs3GmNOA/5ure29u+9MSUmxhYWFQdUlrcBauP9+OOUUSE52YwUFbh9ZWpq7szKUls2GeZMh4N82Fu+DUTOg37jQ/pZICOQWlTElpxh/oHrrmC8+jumZfRk7UAsKIrHCGLPEWpsSzHcENUNmjIkH5gDP7BjGAKy1G6y1G2ufzwfijTH7BfObEiECAbj8cpg8GR5+eNt4aipMmRL6MAaQP61+GAP3On9a6H9LJASy80rqhTEAf6Ca7LwSjyoSkUjV4k39xhgDPAassNbe08g1BwLfWWutMWYQLgCub+lvSoT48Uc46yxYuBD++Ee45ZbW+d2K0uaNi3hsTbm/WeMi0nYFc5flEOBCoNgYs7R27I/AQQDW2oeAs4ArjDFbAD9wrg12jVS8tXq167b/5Zfw5JNw0UWt99sJSVCxuuFxkQjUPdFHWQPhq3uiz4NqRCSStTiQWWvfBsxurrkfuL+lvyERqGtXt19s5kw4/vjW/e30qQ3vIUtXezuJTFkZyQ3uIcvKSPawKhGJRDo6SZrm3/+GjAxITIS5c72poW7jfv40t0yZkOTCmDb0S4Sq27ifnVfCmnI/3RN9ZGUka0O/iOxEgUx2rbrabdLPzoYbb9z1frFls8MflvqNUwATEZGYo0Amjdu40TV7fekluOIKdz5lY3ZsSVGx2r0GBShps3Zse1FW7mdKTjGAZslEpB4dLi4NKy11e8TmzYMZM+CBB2CPXeR3taQQ2YnaXohIU2mGTBrWrh38/LM7n3LEiN1fr5YUIjtR2wsRaSrNkEl9ixe7fWPdu0NxcdPCGDTeekItKaQNa6y9hdpeiMiOFMjEsRZuvx1OOMEdhwTuTMqmSp/qWlBsTy0ppI3LykjGF1///4/U9kJEGqIlS3FLk7/5DTz9NJx3Hlx2WfO/Qy0pRHaithci0lRBHy4eDjpcvBWtWweZmfD22zBtmmttYXbZ71dERES2E4rDxTVDFglao39XY/73P/jkE9f4dZxms0RCLbeoTDNkIrJbCmRe86p/18qV0Ls3pKS4cyn32it8vyXSRqkPmYg0lTb1e82L/l0PPABHHAEvvOBeK4yJhIX6kIlIUymQea01+3e99RakpsLVV7t2FhkZof8NEdlKfchEpKkUyLzWWv27FiyAtDT4739dx/0//AG6dAntb4hIPepDJiJNpUDmtdbq3/Xkk1BT455b62bLRCSs1IdMRJpKgcxr/cbBqBmQ0BMw7nHUjNBt6C8vd49XXQUdO7pmr+3bu9kyEQmrsQN7MD2zLz0SfRigR6KP6Zl9taFfRHaiPmSx7Jln4MorIS8PjjsOCgpg0SIXxlJTva5OREQkJqgPmTSspgZuugluvdWFr9693XhqqoKYiIhIBFIgizWbNsHEifD883DJJfDgg26JUkRERCKW9pDFmpkzXX+x7Gz3XGFMREQk4mmGLFZs2eLaWVx1leu+P3iw1xWJiIhIE2mGLBbMmwdHHgmrV7u7KBXGREREoooCWTSzFu6+G8aMgYQEF8ZEREQk6mjJMlpVVbnlyZkz4ayzXOPXTp0av37ZbHc+ZkWpOwUgfWp4Dy8XEcAdMJ6dV8Kacj/dE31kZSSrD5mI7ESBLFrdeqsLY3/6//buPubu8Y7j+Puj2qqhrYd5aDvKjNgY3R1RE+qpRULFzCqz1WbMxh4sZKSJLBIxk0wiW7JJmYctmOqqhli1lSVTrZuhHlbaklCmphSZ0dZ3f/yum5/T097n9Jz7XOfc5/NKTu7f77qu3znX91znOvf3/B7OmQFXXAFbbWZn51N/hnt+/MmPmK99RqXvQQAACn1JREFUuVgHJ2VmA2jOP1dx2eylH//A+Kq33+ey2UsBnJSZ2af4kGWnuvhiuOuuIjHbXDIGxZ6xdRU/Zrzu/aLczAbMNQ8s+zgZ6/P+ug1c88CyTD0ys3blhKyTLFwIU6YU3zW2ww5w2mm1bbf2lfrKzawpXn37/brKzax7OSHrFDNnwuTJxZWUa9bUt+3IsfWVm1lT7DFqRF3lZta9nJC1uw0bisOT554LxxxT/B7l2DoTqWMvh6EV/wCGjijKzWzAXDJlP0YM/fTVzyOGDuGSKftl6pGZtSsnZO3u4ouLr7a44AK4997i6y3qddAZcPJ1MHIcoOLvydf5hH6zAXbqIWO46rQDGTNqBALGjBrBVacd6BP6zWwjiojcfdhIT09P9Pb25u5Ge1ixAubNg/PPz90TMzMzq0LSYxHR08h9eA9ZO3r0UbjoouKLX/fZx8mYmZnZIOeErN3ceScceSTMmQOrV+fujZmZmbWAE7J28fDDcPzxcMYZMGECLFkCu+6au1dmZmbWAv6m/nawaBEcdRSsX1/8HuWVV8Iuu+TulZmZmbWI95C1g4cego8++mR90aJsXTEzM7PWc0LWDiZNguHDi71jw4YV62ZmZtY1fMiyHUycCPPnF3vKJk0q1s3MzKxrOCFrFxMnOhEzMzPrUj5kaWZmZpaZEzIzMzOzzJyQmZmZmWXWUEIm6QRJyyQtl3Rplfrhku5I9Ysl7dXI45mZmZkNRluckEkaAvwWOBE4ADhT0gEVzc4B3oqIzwPXAldv6eOZmZmZDVaN7CE7FFgeESsj4kPgdmBqRZupwM1peRZwrCQ18JhmZmZmg04jCdkY4OXS+iuprGqbiFgPrAV2auAxzczMzAadtjmpX9J5knol9b7xxhu5u2NmZmbWMo0kZKuAcaX1samsahtJWwMjgTer3VlEXB8RPRHRs4t/WNvMzMy6SCMJ2aPAvpLGSxoGTAPmVrSZC0xPy6cDCyIiGnhMMzMzs0Fni386KSLWS7oQeAAYAtwYEc9IugLojYi5wA3ArZKWA2sokjYzMzMzK2notywj4j7gvoqyy0vL/wO+3shjmJmZmQ12bXNSv5mZmVm3ckJmZmZmlpkTMjMzM7PMnJCZmZmZZeaEzMzMzCwzJ2RmZmZmmTkhMzMzM8vMCZmZmZlZZk7IzMzMzDJzQmZmZmaWmRMyMzMzs8yckJmZmZll5oTMzMzMLDMnZGZmZmaZOSEzMzMzy8wJmZmZmVlmiojcfdiIpHeBZbn7kcHOwH9ydyIDx91dHHd3cdzdpVvj3i8itm/kDrZuVk+abFlE9OTuRKtJ6nXc3cNxdxfH3V0cd3eR1NvoffiQpZmZmVlmTsjMzMzMMmvXhOz63B3IxHF3F8fdXRx3d3Hc3aXhuNvypH4zMzOzbtKue8jMzMzMukaWhEzSjpLmSXoh/R1dpc3BkhZJekbSU5K+Uaq7SdKLkp5It4NbG0F9JJ0gaZmk5ZIurVI/XNIdqX6xpL1KdZel8mWSprSy342qIe6fSXo2je98SXuW6jaUxndua3vemBriPlvSG6X4vleqm57mxQuSpre2542pIe5rSzE/L+ntUl0nj/eNklZLenoT9ZJ0XXpenpI0oVTXyePdX9zfTPEulfSwpC+X6l5K5U804+q0Vqoh7kmS1pZez5eX6jY7R9pZDXFfUor56TSnd0x1HTneksZJWpj+Tz0j6SdV2jRvfkdEy2/Ar4BL0/KlwNVV2nwB2Dct7wG8BoxK6zcBp+fo+xbEOgRYAewNDAOeBA6oaPND4HdpeRpwR1o+ILUfDoxP9zMkd0xNjPtoYNu0/IO+uNP6e7ljGMC4zwZ+U2XbHYGV6e/otDw6d0zNirui/Y+AGzt9vFPfjwQmAE9vov4k4H5AwGHA4k4f7xrjPrwvHuDEvrjT+kvAzrljGKC4JwF/rVJe1xxpt1t/cVe0PRlY0OnjDewOTEjL2wPPV3k/b9r8znXIcipwc1q+GTi1skFEPB8RL6TlV4HVwC4t62HzHAosj4iVEfEhcDtF/GXl52MWcKwkpfLbI+KDiHgRWJ7urxP0G3dELIyI/6bVR4CxLe7jQKhlvDdlCjAvItZExFvAPOCEAepns9Ub95nAbS3p2QCLiL8DazbTZCpwSxQeAUZJ2p3OHu9+446Ih1NcMHjmdy3jvSmNvDdkV2fcg2J+R8RrEfF4Wn4XeA4YU9GsafM7V0K2a0S8lpb/Dey6ucaSDqX4RLGiVHxl2j14raThA9TPZhgDvFxaf4WNB/TjNhGxHlgL7FTjtu2q3r6fQ/Epo882knolPSJpo4S9jdUa99fS63eWpHF1btuOau57OjQ9HlhQKu7U8a7Fpp6bTh7velXO7wD+JukxSedl6tNAmijpSUn3S/piKuuK8Za0LUXicVepuOPHW8WpRIcAiyuqmja/B+yb+iU9COxWpWpGeSUiQtImL/VMmeatwPSI+CgVX0aRyA2juNT058AVzei3tZ6ks4Ae4KhS8Z4RsUrS3sACSUsjYkX1e+g49wC3RcQHkr5PsXf0mMx9aqVpwKyI2FAqG8zj3dUkHU2RkB1RKj4ijfdngXmS/pX2wAwGj1O8nt+TdBIwB9g3c59a6WTgHxFR3pvW0eMtaTuKBPOnEfHOQD3OgO0hi4jjIuJLVW53A6+nRKsv4Vpd7T4k7QDcC8xIuwL77vu1tHvwA+APtPdhvFXAuNL62FRWtY2krYGRwJs1btuuauq7pOMokvRT0ngCEBGr0t+VwEMUn0w6Qb9xR8SbpVhnAl+pdds2Vk/fp1FxOKODx7sWm3puOnm8ayLpIIrX+NSIeLOvvDTeq4G/0N7v4XWJiHci4r20fB8wVNLOdMF4J5ub3x033pKGUiRjf4qI2VWaNG1+5zpkORfou+JgOnB3ZQNJwygG7paImFVR15fMieL8s6pXfbSJR4F9JY1PMU2jiL+s/HycTnEyZKTyaSquwhxP8SlrSYv63ah+45Z0CPB7imRsdal8dN9h6PRG9lXg2Zb1vDG1xL17afUUivMSAB4AJqf4RwOTU1knqOV1jqT9KU5wXVQq6+TxrsVc4NvpaqzDgLXplI1OHu9+SfocMBv4VkQ8Xyr/jKTt+5Yp4m7n9/C6SNot/W/qO91mK4oP2DXNkU4maSTFkY67S2UdO95pHG8AnouIX2+iWfPm95ZcedDojeL8qPnAC8CDwI6pvAeYmZbPAtYBT5RuB6e6BcBSikH9I7BdjjjqiPckiqszVlDs7YPiEOspaXkb4E6Kk/aXAHuXtp2RtlsGnJg7libH/SDweml856byw9P4Ppn+npM7libHfRXwTIpvIbB/advvptfBcuA7uWNpZtxp/RfALyu26/Txvo3iKvB1FOeJnAOcD5yf6gX8Nj0vS4GeQTLe/cU9E3irNL97U/neaayfTPNgRu5Ymhz3haX5/QhweGnbjeZIp9z6izu1OZviQrTydh073hSH2QN4qvQ6Pmmg5re/qd/MzMwsM39Tv5mZmVlmTsjMzMzMMnNCZmZmZpaZEzIzMzOzzJyQmZmZmWXmhMzMzMwsMydkZmZmZpk5ITMzMzPL7P87mqtRCuVlvAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,8))\n",
    "# Plot the graph\n",
    "ax.plot(X,y,'o',label='data')\n",
    "ax.plot(x,y_fitted, 'r--.', label='OLS')\n",
    "# Set the legend\n",
    "ax.legend(loc='best')\n",
    "ax.axis((-0.25,2,-1,20)) ## Set the area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 高次模型的回归\n",
    "\n",
    "$\\Large Y = \\beta_0 + \\beta_1 X + \\beta_2X^2 + ... + \\beta_nX^n$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample:  \n",
    "$\\large Y = 1 + 0.1 \\, X + 10 \\, X^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample2 = 100\n",
    "x2 = np.linspace(0,10,nsample2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.10101</td>\n",
       "      <td>0.010203</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.20202</td>\n",
       "      <td>0.040812</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.30303</td>\n",
       "      <td>0.091827</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.40404</td>\n",
       "      <td>0.163249</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     0        1         2\n",
       "0  1.0  0.00000  0.000000\n",
       "1  1.0  0.10101  0.010203\n",
       "2  1.0  0.20202  0.040812\n",
       "3  1.0  0.30303  0.091827\n",
       "4  1.0  0.40404  0.163249"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X2 = np.column_stack((x2, x2**2))\n",
    "X2 = sm.add_constant(X2)\n",
    "pd.DataFrame(X2).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.   0.1 10. ]\n"
     ]
    }
   ],
   "source": [
    "# Set the beta_0, beta_1, beta_2 as 1, 0.1, 10\n",
    "beta2 = np.array([1,0.1,10])\n",
    "print(beta2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "y2 = np.dot(X,beta) + e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>x</th>\n",
       "      <th>kx</th>\n",
       "      <th>error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.667889</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>-1.437533</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.069557</td>\n",
       "      <td>0.10101</td>\n",
       "      <td>1.112131</td>\n",
       "      <td>0.452448</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.426198</td>\n",
       "      <td>0.20202</td>\n",
       "      <td>1.428324</td>\n",
       "      <td>1.189455</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.278692</td>\n",
       "      <td>0.30303</td>\n",
       "      <td>1.948577</td>\n",
       "      <td>-0.755288</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.481142</td>\n",
       "      <td>0.40404</td>\n",
       "      <td>2.672891</td>\n",
       "      <td>2.176870</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          y        x        kx     error\n",
       "0  1.667889  0.00000  1.000000 -1.437533\n",
       "1  0.069557  0.10101  1.112131  0.452448\n",
       "2  0.426198  0.20202  1.428324  1.189455\n",
       "3  0.278692  0.30303  1.948577 -0.755288\n",
       "4  1.481142  0.40404  2.672891  2.176870"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e2 = np.random.normal(size=nsample2)\n",
    "pd.DataFrame({\n",
    "    \"y\": y2,\n",
    "    \"x\": x2,\n",
    "    \"kx\": np.dot(X2,beta),\n",
    "    \"error\": e2\n",
    "}).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "model2 = sm.OLS(y2,X2)\n",
    "results2 = model2.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.65396289  0.12220614 10.00800226]\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       1.000\n",
      "Model:                            OLS   Adj. R-squared:                  1.000\n",
      "Method:                 Least Squares   F-statistic:                 4.256e+06\n",
      "Date:                Sun, 26 Apr 2020   Prob (F-statistic):          1.78e-240\n",
      "Time:                        22:41:49   Log-Likelihood:                -143.79\n",
      "No. Observations:                 100   AIC:                             293.6\n",
      "Df Residuals:                      97   BIC:                             301.4\n",
      "Df Model:                           2                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.6540      0.304      2.149      0.034       0.050       1.258\n",
      "x1             0.1222      0.141      0.869      0.387      -0.157       0.401\n",
      "x2            10.0080      0.014    735.335      0.000       9.981      10.035\n",
      "==============================================================================\n",
      "Omnibus:                        1.892   Durbin-Watson:                   1.915\n",
      "Prob(Omnibus):                  0.388   Jarque-Bera (JB):                1.373\n",
      "Skew:                          -0.262   Prob(JB):                        0.503\n",
      "Kurtosis:                       3.236   Cond. No.                         144.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "print(results2.params)\n",
    "print(results2.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Find the fittedvalues get the y\n",
    "y_fitted_new = results2.fittedvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 2, -1, 50)"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAHWCAYAAABAA0zqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde3RU1d3/8fchBBMQiVqUmwregraA0KAiVlFaoqKIaFHbem1r76ZVU0WtpvBYbdMb/dl6q32srU8VLUYRNbVRWh/hUS5BwGq8UK0EqBYNcplICOf3xyFcNEBCZubMJO/XWqyZ2TmZ88XVpZ/u/T17B2EYIkmSpNToFHcBkiRJ7ZlhS5IkKYUMW5IkSSlk2JIkSUohw5YkSVIKGbYkSZJSqHNLLgqC4E1gDdAIbAzDsCgIgn2AB4D+wJvAxDAM309NmZIkSdmpNTNbJ4VheFQYhkWbP18DVIVheBhQtfmzJEmSttGWZcQzgd9vfv97YHzby5EkSWpfWhq2QuAvQRDMD4Lgss1j+4dhuGLz+5XA/kmvTpIkKcu1qGcLOD4Mw9ogCPYDngqC4JVtfxiGYRgEQbPn/mwOZ5cBdOvW7dMDBw5sU8GSJEnpMH/+/P+EYdizrd/TorAVhmHt5td3giB4GDga+HcQBL3DMFwRBEFv4J0d/O6dwJ0ARUVF4bx589pasyRJUsoFQfBWMr5nl8uIQRB0C4Kge9N7YAywBHgUuGjzZRcBjySjIElqqYrqWkbe8jQDrpnJyFuepqK6Nu6SJOljWjKztT/wcBAETdf/TxiGTwZBMBeYFgTBl4G3gImpK1OStldRXcuk6YtJNDQCUFuXYNL0xQCMH9o3ztIkaTu7DFthGC4FhjQzvgoYnYqiJGlXyitrtgStJomGRsorawxbkjJKSxvkJSmjLK9LtGpcUus1NDSwbNky6uvr4y4lpfLy8ujXrx+5ubkp+X7DlqSs1Kcgn9pmglWfgvwYqpHap2XLltG9e3f69+/P5naidicMQ1atWsWyZcsYMGBASu7h2YiSslJpcSH5uTnbjeXn5lBaXBhTRVL7U19fz7777ttugxZAEATsu+++KZ29c2ZLUlZq6ssqr6xheV2CPgX5lBYX2q8lJVl7DlpNUv13NGxJylrjh/Y1XEkdTFlZGXvuuSdXXXVVsz+vqKjg8MMP58gjj0xzZTvmMqIkSUqKTNj7rqKign/84x9pv+/OGLYkSVKbNe19V1uXIGTr3nfJCFw33XQThx9+OMcffzw1NTUA3HXXXQwfPpwhQ4Zw9tlns379embPns2jjz5KaWkpRx11FG+88Uaz16WbYUuSJLXZzva+a4v58+dz//33s3DhQh5//HHmzp0LwIQJE5g7dy4vvvgiRxxxBHfffTfHHXcc48aNo7y8nIULF3LIIYc0e1262bMlSZLaLFV73z377LOcddZZdO3aFYBx48YBsGTJEq6//nrq6upYu3YtxcXFzf5+S69LJWe2JElSm+1oj7tU7X138cUXc+utt7J48WJuvPHGHW7d0NLrUsmwJUmS2ixVe9+dcMIJVFRUkEgkWLNmDTNmzABgzZo19O7dm4aGBu67774t13fv3p01a9Zs+byj69LJsCVJktps/NC+3DxhEH0L8gmAvgX53DxhUJu3Zxk2bBjnnnsuQ4YM4dRTT2X48OEATJkyhWOOOYaRI0cycODALdefd955lJeXM3ToUN54440dXpdOQRiGabtZUVFROG/evLTdT5Ik7b6XX36ZI444Iu4y0qK5v2sQBPPDMCxq63c7syVJkpRChi1JkqQUMmxJkiSlkGFLkiQphQxbkiRJKWTYkiRJSiHDliRJyhplZWX89Kc/3eHPKyoq+Mc//pHGinbNsCVJkpJj0TT4xaegrCB6XTQt7SUYtiRJUvu0aBrMuBxWvw2E0euMy5MSuG666SYOP/xwjj/+eGpqagC46667GD58OEOGDOHss89m/fr1zJ49m0cffZTS0lKOOuoo3njjjWavSzfDliRJaruqydCQ2H6sIRGNt8H8+fO5//77WbhwIY8//jhz584FYMKECcydO5cXX3yRI444grvvvpvjjjuOcePGUV5ezsKFCznkkEOavS7dOqf9jpIkqf1Zvax14y307LPPctZZZ9G1a1cAxo0bB8CSJUu4/vrrqaurY+3atRQXFzf7+y29LpWc2ZIkSW3Xo1/rxtvo4osv5tZbb2Xx4sXceOON1NfXt+m6VDJsSZKktht9A+Tmbz+Wmx+Nt8EJJ5xARUUFiUSCNWvWMGPGDADWrFlD7969aWho4L777ttyfffu3VmzZs2Wzzu6Lp0MW5Ikqe0GT4QzfgU9DgCC6PWMX0XjbTBs2DDOPfdchgwZwqmnnsrw4cMBmDJlCscccwwjR45k4MCBW64/77zzKC8vZ+jQobzxxhs7vC6dgjAM03azoqKicN68eWm7n6T2raK6lvLKGpbXJehTkE9pcSHjh/aNuyyp3Xj55Zc54ogj4i4jLZr7uwZBMD8Mw6K2frcN8pKyUkV1LZOmLybR0AhAbV2CSdMXAxi4JGUUlxElZaXyypotQatJoqGR8sqamCqSpOYZtiRlpeV1iVaNS1JcDFuSslKfgvxWjUtSXAxbkrJSaXEh+bk5243l5+ZQWlwYU0WS1Dwb5CVlpaYmeJ9GlJTpDFuSstb4oX0NV1IHU1ZWxp577slVV13V7M8rKio4/PDDOfLII9Nc2Y65jChJkpJi5tKZjHloDIN/P5gxD41h5tKZaa+hoqKCf/zjH2m/784YtiRJUpvNXDqTstllrFi3gpCQFetWUDa7LCmB66abbuLwww/n+OOPp6Ym2t7lrrvuYvjw4QwZMoSzzz6b9evXM3v2bB599FFKS0s56qijeOONN5q9Lt0MW5Ikqc2mLphKfeP2hzzXN9YzdcHUNn3v/Pnzuf/++1m4cCGPP/44c+fOBWDChAnMnTuXF198kSOOOIK7776b4447jnHjxlFeXs7ChQs55JBDmr0u3ezZkiRJbbZy3cpWjbfUs88+y1lnnUXXrl0BGDduHABLlizh+uuvp66ujrVr11JcXNzs77f0ulRyZkuSJLVZr269WjXeVhdffDG33norixcv5sYbb6S+vr5N16WSYUuSJLVZybAS8nLythvLy8mjZFhJm773hBNOoKKigkQiwZo1a5gxYwYAa9asoXfv3jQ0NHDfffdtub579+6sWbNmy+cdXZdOLiNKkqQ2G3vwWCDq3Vq5biW9uvWiZFjJlvHdNWzYMM4991yGDBnCfvvtx/DhwwGYMmUKxxxzDD179uSYY47ZErDOO+88vvrVr/KrX/2Khx56aIfXpVMQhmHablZUVBTOmzcvbfeTJEm77+WXX+aII46Iu4y0aO7vGgTB/DAMi9r63S4jSpIkpZBhS5IkKYUMW5IkSSlk2JIkSUohw5YkSVIKGbYkSZJSyLAlSZIy2rJlyzjzzDM57LDDOOSQQygpKWHDhg3MmjWL008//WPXP/bYYwwdOpQhQ4Zw5JFHcscdd8RQ9VaGLUmSlLHCMGTChAmMHz+e1157jVdffZW1a9dy3XXXNXt9Q0MDl112GTNmzODFF1+kurqaUaNGpbfoj3AHeUmSlDxz5sCsWTBqFIwY0eave/rpp8nLy+OSSy4BICcnh1/84hcMGDCAk0466WPXr1mzho0bN7LvvvsCsMcee1BYWNjmOtrCsCVJklqmuRmiiRPhm9+E9eth5EhYtAg2bYJOnWDwYCgpgYsvhv/8B845Z/vfnTVrl7d86aWX+PSnP73d2F577cWBBx7I66+//rHr99lnH8aNG8dBBx3E6NGjOf300zn//PPp1Cm+xTyXESVJUnKsXh0FLYheV6+OpYzf/va3VFVVcfTRR/PTn/6USy+9NJY6mjizJUmSWmZnM1Fdu8J998Ho0bBhA3TpEn1uWkr8xCdaNJP1UUceeSQPPfTQdmMffPAB//rXvzj00EP5y1/+0uzvDRo0iEGDBnHBBRcwYMAA7rnnnlbfO1mc2ZIkSckxYgRUVcGUKdFrEnq2Ro8ezfr167n33nsBaGxs5Morr+Tiiy+ma9euH7t+7dq1zNom1C1cuJCDDjqozXW0hWFLkiQlz4gRMGlSUoIWQBAEPPzwwzz44IMcdthhHH744eTl5fGjH/0IgKqqKvr167flT3V1NT/5yU8oLCzkqKOO4sYbb4x1VgtcRpQkSRnugAMOYMaMGR8bHzVqFIlE4mPjn/nMZ9JRVosZtiRlrYrqWsora1hel6BPQT6lxYWMH9o37rIkaTuGLUlZqaK6lknTF5NoaASgti7BpOmLAQxckjKKPVuSslJ5Zc2WoNUk0dBIeWVNTBVJUvMMW5Ky0vK6j/dp7Gxc0u4JwzDuElIu1X9Hw5akrNSnIL9V45JaLy8vj1WrVrXrwBWGIatWrSIvLy9l97BnS1JWKi0u3K5nCyA/N4fS4njPQJPak379+rFs2TLefffduEtJqby8PPr165ey7zdsScpKTU3wPo0opU5ubi4DBgyIu4ysZ9iSlLXGD+1ruJKU8ezZkiRJSiHDliRJUgoZtiRJklLIsCVJkpRChi1JkqSPqq9P2lcZtiRJkrYVhvCVryTt6wxbkiRJ25oxA+67L2lfZ9iSJEna1hlnwMMPJ+3rDFuSJEkAL70Er70GQQDjxyfta91BXpIk6d134fTToVs3WLQIOiVvPqrF3xQEQU4QBNVBEDy2+fOAIAieD4Lg9SAIHgiCoEvSqpIkSUqXDRvg7LNh5Ur43e+SGrSgdcuIJcDL23z+MfCLMAwPBd4HvpzMwiRJklIuDOGb34Rnn42C1tFHJ/0WLQpbQRD0A8YCv938OQBOBh7afMnvgeQtbkqSJKXD//wP3H03XHcdnH9+Sm7R0p6tXwLfB7pv/rwvUBeG4cbNn5cBfZNcmyRJUmpNmAC/+AVcfnnKbrHLma0gCE4H3gnDcP7u3CAIgsuCIJgXBMG8d999d3e+QpIkKbn++U9YvRry8+G73016n9a2WvLNI4FxQRC8CdxPtHw4FSgIgqBpZqwfUNvcL4dheGcYhkVhGBb17NkzCSVLkiS1wfvvQ3ExnHlm1LOVYrsMW2EYTgrDsF8Yhv2B84CnwzD8IvAMcM7myy4CHklZlZIkScmwcSOcey68+SbcdFO0p1aKtWXO7GrgiiAIXifq4bo7OSVJkiSlyJVXwlNPwe23w8iRabllqzY1DcNwFjBr8/ulQPKfj5QkSUqFe+6BX/0Kvvc9uPTStN3W43okSVLH8NnPwhVXwE9+ktbbGrYkSVL79u670NgI/frBz34GndN7WqFhS5IktV8ffAAnnQSXXBJbCYYtSZLUPjU2whe/CK+8AhddFFsZ6Z1HkyRJSpfrroPHHoNf/xpGj46tDGe2JElS+/OHP8CPfwzf+EZ00HSMDFuSJKn9OewwOO88mDo17kpcRpQkSe1IfT3k5cGxx0Z/MoAzW5IkqX1Yty7aFX7KlLgr2Y4zW5KyVkV1LeWVNSyvS9CnIJ/S4kLGD+0bd1mS0mHRNKiaDKuXQY9+cNL1cNOfobrasCVJyVBRXcuk6YtJNDQCUFuXYNL0xQAGLqm9WzQNZlwODYno8+q34cqvwNProLwcTjst3vo+wmVESVmpvLJmS9BqkmhopLyyJqaKJKVN1eStQQvgpYYoaA0viA6azjDObEnKSsvrEq0al9SOrF4Wvb69Ed5shI0h9M+BMZsgCOKtrRmGLUlZqU9BPrXNBKs+BfkxVCMprXr0gyX/hHvXQyOQA1zQFfY9IO7KmuUyoqSsVFpcSH5uznZj+bk5lBYXxlSRpLQZfQO8GcBGICQKXMs6ReMZyLAlKSuNH9qXmycMom9BPgHQtyCfmycMsjle6gg+eTa81y96HwCdA7jgahg8MdaydsRlRElZa/zQvoYrqaMJQygpgYU1cMUV8IlPwKhRMGJE3JXtkGFLkiRlj5//PDpY+qqrom0esoDLiJIkKXu88w58/vPRIdNZwpktSZKU+TZtgk6dopDV2Bi9zxLZU6kkSeqYXnsNjjoKFi6MPufk7Pz6DOPMliRJylzvvgunngqrV0P37nFXs1sMW5IkKTMlEjBuHNTWwjPPwCGHxF3RbjFsSZKkzNPYCF/6Ejz/PDz0EBx7bNwV7TZ7tiRJUub58ENYvx5+9jOYMCHuatrEmS1JkpRZNm2Crl3hscey6qnDHcn+v4EkSWo/HnkERo6MGuNzciAI4q6ozQxbkiQpM7zwApx/fjSz1a1b3NUkjWFLkiTFb+lSOOMM6NULZsyIlhHbCcOWJEmK13vvwWmnQUMDPPEE7Ldf3BUllWFLkiTF64MPoEuXqF+rsDDuapLOpxElSVI8Nm2KGuD794fq6qw7hqelnNmSJEnxuO46uOAC2Lix3QYtMGxJkqQ43HEH3HIL7Llnuw5aYNiSJEnp9vjj8M1vRk3xt97aLvbS2hnDliRJSp8FC2DiRBgyBB54ADq3//Zxw5YkSUqfDz6Agw+OjuLZc8+4q0mL9h8nJUlS/Bobo96sUaNg4cJ2ceZhS3Wcv6kkSYrHhg1wyinwi19EnztQ0ALDliRJSqUwhK98Bf76V9h337iriYVhS5IkpU5ZGfzhD/DDH8KFF8ZdTSzs2ZKUtSqqaymvrGF5XYI+BfmUFhcyfmjfuMuSOqZF06BqMqxeBj36wegbYP46mDwZLrkEfvCDuCuMjWFLUlaqqK5l0vTFJBoaAaitSzBp+mIAA5eUboumwYzLoSERfV79dvR53dioV+uOO9r9Xlo74zKipKxUXlmzJWg1STQ0Ul5ZE1NFUgdWNXlr0ALYFEafez4PM2dCbm58tWUAw5akrLS8LtGqcUkptHrZ1vcfbILb1sEbG6PxDvbkYXP8JyApK/UpyG/VuKQU6tEvel26Ee5aB3WboFuwdbyDM2xJykqlxYXk525/eG1+bg6lxYUxVSR1YKNvgGU58Mf1sDaEEAi7ROMybEnKTuOH9uXmCYPoW5BPAPQtyOfmCYNsjpfi8Mmz4W97RiELYBOQVwyDJ8ZZVcbwaURJWWv80L6GKykThCEMHAJvrYRNm6BLFzjv8riryhiGLUmStPvWro0OlH70Ufi//4NZs6LzD0eMiLuyjOEyoiRJ2j2//CUMHgzLl0f7aI0YAZMmGbQ+wrAlSZJa79574Xvfg6FDYb/94q4moxm2JElS68yYAZdeCiefDPfdB53tStoZw5YkSWq5OXNg4sRoRquiAvLy4q4o4xm2JElSyxUWwuc/D48/Dt27x11NVnDeT5Ik7drbb0PPnrDPPlG/llrMmS1JkrRzy5fDCSfARRfFXUlWMmxJkqQde/99KC6G//wHrroq7mqyksuIkiSpeevWwemnw6uvRj1aw4fHXVFWMmxJkqTmffWr0a7wDz4Io0fHXU3WMmxJkqTmXXstnHYaTJgQdyVZzZ4tSZK0VRhCZWX0+qlPwZe+FHdFWc+wJUmStpoyBU45JdolXklh2JIkSZFf/xpuvBEuvhjOOCPuatoNw5YkSYI//Qm+8x0YNw7uuguCIO6K2g3DliRJHd2KFdHB0iecAPff78HSSeY/TUmSOrrevaNDpY89FvLz466m3XFmS5KkjmrJEpg5M3pfXAw9esRbTzvlzJYkSR3RP/8JY8ZAbi7U1EBeXtwVtVuGLUmSOpp//xs+9zmor4e//MWglWKGLUmSOpK6umjJcMUKqKqKNi5VShm2JEnqSO69F/7xj2jT0mOPjbuaDsGwJUlSR/Kd78BJJ8GgQXFX0mEYtiRlrYrqWsora1hel6BPQT6lxYWMH9o37rKkzLBoGlRNhtXLoHtfWHIIXP0TGDjQoJVmhi1JWamiupZJ0xeTaGgEoLYuwaTpiwEMXNKiaTDjcmhIRAdKP/g6PP8P2Hsv+Pkf466uw3GfLUlZqbyyZkvQapJoaKS8siamiqQMUjU5Clpvb4Q/rofnN8DRXaBfddyVdUjObEnKSsvrEq0alzqU1cuioHXPetgEBMAnO8MHtXFX1iHtcmYrCIK8IAheCILgxSAIXgqC4IebxwcEQfB8EASvB0HwQBAEXVJfriRF+hQ0f6TIjsalDqVHP3hzYxS0mrzVGI0r7VqyjPghcHIYhkOAo4BTgiA4Fvgx8IswDA8F3ge+nLoyJWl7pcWF5OfmbDeWn5tDaXFhTBVJGeTkH8Ch3aL1qwDIAQ7tCqNviLmwjmmXYSuMrN38MXfznxA4GXho8/jvgfEpqVCSmjF+aF9unjCIvgX5BEDfgnxunjDI5njpwQeh5Da46OfwzQFw0h7R67duh8ET466uQ2pRz1YQBDnAfOBQ4NfAG0BdGIYbN1+yDGj233BBEFwGXAZw4IEHtrVeSdpi/NC+hitpWxUV8IUvRJuVDj4Hjrs47opEC59GDMOwMQzDo4B+wNHAwJbeIAzDO8MwLArDsKhnz567WaYkSdqpxx6DiROhqAgefxz23DPuirRZq7Z+CMOwDngGGAEUBEHQNDPWD/ARB0mS4vDUU3D22TBkCDzxBHTvHndF2kZLnkbsGQRBweb3+cDngJeJQtc5my+7CHgkVUVKkqSdOOggGDMGKiuhoCDuavQRLZnZ6g08EwTBImAu8FQYho8BVwNXBEHwOrAvcHfqypQkSR/z2mvRDvGHHx4dLL3PPnFXpGbsskE+DMNFwNBmxpcS9W9JkqR0e+45KC6Ga6+N/ihjeVyPJEnZ5vnn4dRToU8fuOSSuKvRLhi2JEnKJvPnRzNaPXvC009D795xV6RdMGxJkpQt1q+HsWOjJvinn4Z+Hr+TDTyIWpKkbNG1K9xzT9QQf9BBcVejFnJmS5KkTPfKK/DnP0fvTzkFDj443nrUKs5sSZKUyV57DU4+GTp1ioJWt25xV6RWMmxJkpSp3ngDTjoJNm6EWbMMWlnKsCVJUiZ6661oRiuRgGeegSOPjLsi7SbDliRJmejBB+GDD6CqCgYPjrsatYEN8pIkZZIwjF6vvBKWLIFhw+KtR21m2JIkKVP8+98wahQsWgRBAH37xl2RksBlREmSMsG778Lo0fDPf8Lq1XFXoyQybEmSFLf33oPPfS56+nDmTPjMZ+KuSElk2JIkKU51dTBmDLz8MsyYET2BqHbFni1JkuLUpUt0mPTDD0ehS+2OM1uSslZFdS3llTUsr0vQpyCf0uJCxg+1oVgZatE0qJoMq5dBj34w4vtwxJnQowc8+mjUEK92yZktSVmporqWSdMXU1uXIARq6xJMmr6YiurauEuTPm7RNJhxOax+Gwjh3X/B+V+BE4+GTZsMWu2cYUtSViqvrCHR0LjdWKKhkfLKmpgqknaiajI0JKL3DSHcvx7eaoDB66MzD9WuuYwoKSstr0u0alyK1epl0eubG2FGAt4LYXweHLwm3rqUFsZpSVmpT0F+q8alWPXoB29vhHvXR0GrE7BPp2hc7Z5hS1JWKi0uJD83Z7ux/NwcSosLY6pI2onRN8Db2/wnNyT6PPqG2EpS+hi2JGWl8UP7cvOEQfQtyCcA+hbkc/OEQT6NqMyzZg089TZ8sRQ6BxAQvV5wNQyeGHd1SoMgbDrwMg2KiorCefPmpe1+kiTFqq4OTjkF5s2DOXNg40aYNSs6/3DEiLir0y4EQTA/DMOitn6PDfKSJKXCf/4TbVK6ZAk8+CAMHx6NG7I6HMOWJEnJtnIlfPaz0VmHjzwCp54ad0WKkWFLkqRke+UVWLEiOlTasw47PMOWJEnJsn49dO0a9WS9+SZ07x53RcoAPo0oSVIy1NTAwIFw333RZ4OWNnNmS5Kktlq8OOrRCkMYNCjuapRhnNmSJKktFiyIlg07d4a//x0GD467ImUYw5YkSbtr+fKoAb579yhoDRwYd0XKQIYtSZJ2V58+8F//Bc8+C4ccEnc1ylD2bEmS1FqVlfCJT8CnPw3f/nbc1SjDObMlSVJrPPIIjBsHV18ddyXKEoYtSZJa6oEH4Oyz4aijoiN4pBYwbEmS1BL33ANf+AIcdxw89RTsvXfcFSlLGLYkSdqVMISHH4bRo+HJJ2GvveKuSFnEBnlJknamvh7y8qIlRIjeS63gzJYkSTty001w7LGwenUUsgxa2g2GLUmSPioM4brr4Prrox3hu3WLuyJlMZcRJUnaVhjCFVfAL38Jl10Gt90GnZyb0O7zfz2SJG1rypQoaF1+Odx+u0FLbebMliRJ27rkEsjPh6uugiCIuxq1A8Z1SZIaGuDXv4bGRjjgACgtNWgpaZzZkpS1KqprKa+sYXldgj4F+ZQWFzJ+aN+4y1K2WDQNqibDqrehIoSX1sChh0JxcdyVqZ0xbEnKShXVtUyavphEQyMAtXUJJk1fDGDg0q4tmgYzLodX18BjCXgvhDO6Q+/VcVemdshlRElZqbyyZkvQapJoaKS8siamipRVqiZDzRr4w/ooaHUCejZG41KSGbYkZaXldYlWjUvbWb0MXm6AcPPnEHizMRqXksywJSkr9SnIb9W4tMX770OPfjA4N2qmCYAcoH9ONC4lmWFLUlYqLS4kPzdnu7H83BxKiwtjqkhZYe5cKCyEtSfAwd3hwq5w0h7R68HdYfQNcVeodsgGeUlZqakJ3qcR1WJPPglnnw377QfnfhfqR0U9Wgcsi2a0Rt8AgyfGXaXaoSAMw11flSRFRUXhvHnz0nY/SZIAuPde+PKX4VOfgieegF694q5IWSAIgvlhGBa19XtcRpQktW+LF8NFF8GJJ8Lf/mbQUtq5jChJat8GDYKHH4ZTT4U99oi7GnVAzmxJktqfDz+Mzjh87rno8/jxBi3FxrAlSWpfPvgATjsN7rkHFiyIuxrJZURJUjuycmW0XLhkSdQUf8EFcVckGbYkSe3EihUwciS88w489pgHSitjuIwoSWof9tsPxoyBZ54xaCmjOLMlScpuTz0FAwfCAQfA7bfHXY30Mc5sSZKy1733Rs3wV18ddyXSDhm2JEnZJwzhxz/eulmpM1rKYIYtSVJ22bQJvvc9uOYaOP98ePxx2GuvuKuSdsiwJUnKLuvWRU3w3/0u/PGP0KVL3BVJO2WDvCQpO3zwAeTmQvfu0c7w3bpBEKoaQlIAACAASURBVMRdlbRLzmxJkjLfihVwwglw4YXR5z33NGgpaxi2JEmZ7dVX4bjj4PXX4StfibsaqdVcRpQkZa4XXoCxY6NZrGeegeHD465IajXDliQpM23YABMnRj1alZVw2GFxVyTtFsOWJCkzdekC06dDnz7Qq1fc1Ui7zbAlScocYQg/+QmsXw8//CEMGxZ3RVKb2SAvScoMzz0HI0dGm5W+9lq0eanUDjizJUmKX1UVjBkTBazOneFb34JOzgeoffB/yZKyVkV1LSNveZoB18xk5C1PU1FdG3dJ2h2NjfDFiVtnsho3wrRb461JHdrMpTMZ89AY8vrnfToZ32fYkpSVKqprmTR9MbV1CUKgti7BpOmLDVzZ6KU/wxH1kAMERK/1lbBoWsyFqSOauXQmZbPLWLFuRdK+02VESVmpvLKGREPjdmOJhkbKK2sYP7RvTFWpVe6/H3JyYNlkOLEzHNwV3myE/jnQuxGqJsPgiXFXqQ5m6oKp1DfWJ/U7dzmzFQTBAUEQPBMEwT+CIHgpCIKSzeP7BEHwVBAEr21+3TuplUnSTiyvS7RqXBkkDOHGG+H88+G3v4W6t6PxAzrDZ/aIXgFWL4uvRnVYK9etTPp3tmQZcSNwZRiGRwLHAt8KguBI4BqgKgzDw4CqzZ8lKS36FOS3alwZIpGA886DyZPh4ovh0Ueh4IDmr+3RL62lSQC9uiV/T7ddhq0wDFeEYbhg8/s1wMtAX+BM4PebL/s9MD7p1UnSDpQWF5Kfm7PdWH5uDqXFhTFVpF1KJODEE+HBB6O9tH73O9hjDxh9A+R+JCTn5kfjUpqVDCshLycvqd/Zqp6tIAj6A0OB54H9wzBs6h5bCeyf1MokaSea+rLKK2tYXpegT0E+pcWF9mtlsvx8+Oxn4frrYdy4reNNfVlVk6Olwx79oqBlv5ZiMPbgsUDUu/U6ryflO4MwDFt2YRDsCfwNuCkMw+lBENSFYViwzc/fD8PwY31bQRBcBlwGcOCBB376rbfeSkrhkqQs8fDD0L8/DB0adyVSqwRBMD8Mw6K2fk+Ltn4IgiAX+DNwXxiG0zcP/zsIgt6bf94beKe53w3D8M4wDIvCMCzq2bNnW+uVJGWLMISbb4YJE+Cmm+KuRopNS55GDIC7gZfDMPz5Nj96FLho8/uLgEeSX54kKSt9+CFcdBFcey184Qvwxz/GXZEUm5b0bI0ELgAWB0GwcPPYtcAtwLQgCL4MvAW4uC5Jgvffh9NPh9mzYcoUuO46CIK4q5Jis8uwFYbh/xLt6duc0cktR5KU9bp3h549o6cOzzkn7mqk2LmDvCQpOSor4aijYP/9o6Z4Z7MkwLMRJUltFYbw85/DqafCDZv3xjJoSVs4syVJ2n0bNsC3vhUdu3POOfCLX8RdkZRxnNmSJO2eVatgzJgoaF1/PTzwAHTtGndVUsZxZkuStPtWrYL77ou2d5DULMOWJKl15syBYcNg332huho6+58SaWdcRpQktdyvfw2f+Uy0MzwYtKQWMGxJknZt48aoEf7b34bTToMrr4y7IilrGLYkSTtXVxcFrN/8BkpLoz20unePuyopazj/K0lq3pw5MGsW9O8PCxbA3XfDpZfGXZWUdQxbkqSPmzMHTj4ZGhqgSxd49FH47GfjrkrKSoYtSVmrorqW8soaltcl6FOQT2lxIeOH9o27rOwXhnDVN6C+Pvpcn4AZdxu2stTMpTOZumAqK9etpFe3XpQMK2HswWPjLqtDsWdLUlaqqK5l0vTF1NYlCIHaugSTpi+moro27tKy2+rVMPpomP0iBER/coD6Slg0Lebi1Fozl86kbHYZK9atICRkxboVlM0uY+bSmXGX1qEYtiRlpfLKGhINjduNJRoaKa+siamidmDRIigqgr/Ng+I94OKucNIecGFX6N0IVZPjrlCtNHXBVOob67cbq2+sZ+qCqTFV1DG5jCgpKy2vS7RqXC2wbBkkEnBRNzgwJxo7cJv/TKxeFk9d2m0r161s1bhSw5ktSVmpT0F+q8a1Ax9+CE89Fb0/7TR47TUY1L/5a3v0S1tZSo5e3Xq1alypYdiSlJVKiwvJz83Zbiw/N4fS4sKYKspC//oXnHACnHoqLF0ajeXnw+gbIPcjoTV387iySsmwEvJy8rYby8vJo2RYSUwVdUwuI0rKSk1PHfo04m566ik4/3zYsAGmTYODD976s8ETo9eqydHSYY9+UdBqGlfWaHrq0KcR4xWEYZi2mxUVFYXz5s1L2/0kSc245Ra49lo48kj485+h0NlAqTlBEMwPw7Cord/jMqIkdTRBAOedB88/b9CS0sBlREnqCBYuhFWrYPRo+P73o7EgiLcmqYNwZkuS2rt77oERI+CKK2DTpihkGbSktDFsSVJ7VV8Pl10Gl1wCxx0XNcV38l/7Urq5jCgpa3WosxEXTWvd04EffBAdJD1/PkyaBFOmQE7Ojq/Pcp7/p0xm2JKUlZrORmw6sqfpbESg/QWuRdNgxuXQsHl3/NVvR59hx4Gre3c4+mj4wQ/gzDPTU2dMms7/azqWpun8P8DApYzgfLKkrNShzkasmrw1aDVpSHz8rMJNm+BHP4Kamqgn6ze/afdBCzz/T5nPsCUpK3WosxF3dCbhtuPvvQennw7XXQd/+lN66soQnv+nTGfYkpSVOtTZiDs6k7BpfP58GDYMqqrg9tvhxhvTV1sG8Pw/ZTrDlqSs1KHORtzZWYXPPgsjR0ZLiM8+C1/7Wofb1sHz/5TpbJCXlJU61NmIOzursL4evv51uP56+MQn4q0zJp7/p0zn2YiSlG2WLoWrr4bf/hZ69Ii7Gqnd8mxESepo5syBiy6CIUPgr3+Fl1+OuyJJLeAyoiRlg2efjTYp3bgx6smaNg2OPTbuqiS1gDNbkpQNrrkmCloQHbnz2mvx1iOpxQxbkpSpwhASm/cN+/73oUuX6MidLl1g1KhYS5PUci4jSlImqquLnjJ8/3144oloJ/hZs6I/o0bBiBExFyippQxbkpRpnn0WvvQlWL4cfvjDaIYLooBlyJKyjsuIkpQpGhqig6NHjYLcXHjuObj22mjpUFLWMmxJUqZYuxbuuQcuvBCqq+Hoo+OuSFISuIwoSXEKQ3jkETjtNNh77yhkddCd4KX2ypktSYpLXR2cfz6cdRbcfXc0ZtCS2h1ntiQpDts2wd90E1x2WdwVSUoRZ7YkKd1uu80meKkDMWxJUrqNHAmXXmoTvNRBGLYkKdXCEO69F0pKos+DB8Ndd0H37vHWJSktDFuSlEpNTfAXXQQLF249fkdSh2GDvKSsVVFdS3llDcvrEvQpyKe0uJDxQ/um5+aLpkHVZFi9DHr0g9E3wOCJ21/z97/DBRdAbW3UBH/11UnvzZq5dCZTF0xl5bqV9OrWi5JhJYw9eGxS7yGpbQxbkrJSRXUtk6YvJtHQCEBtXYJJ0xcDpD5wLZoGMy6Hhs2zVKvfjj7D1sC1Zg2MHw/77AOzZ6ekN2vm0pmUzS6jvrEegBXrVlA2uwzAwCVlEJcRJWWl8sqaLUGrSaKhkfLKmtTfvGry1qDVpCERjdfWRj1a3bvDY4+ltAl+6oKpW4JWk/rGeqYumJqS+0naPYYtSVlpeV3zvU87Gk+q1cs+PhaG8PelMHAg3H57NHbccSltgl+5bmWrxiXFw7AlKSv1Kchv1XhS9ei3/ef6EP6cgIoEDBsGY9OzhNerW69WjUuKh2FLUlYqLS4kP3f7ZvP83BxKiwtTf/PRN0BuPry9ER5LwK1r4OWN8J3z4Omn4cADU18DUDKshLycvO3G8nLyKBlWkpb7S2oZG+QlZaWmJvhYnkYcPBFefBV+eAM0hNHYjV+DsttTf+9tNDXB+zSilNkMW5Ky1vihfdO31cO2XngBHnwBNnUCGqPtHPY4KP11EAUuw5WU2VxGlKSWWr8erroKRoyAuXOhS5coaHXpEp11KEnNcGZLklrib3+DL38Z3ngDvv51+PGP4aWXYNasKGiNGBF3hZIylGFLknbl3/+G4mLo1w+eeWbrLNaIEYYsSbvkMqIk7Uh1dfS6//4wYwYsWuRyoaRWM2xJ0ketWgUXXhjtmVVZGY197nPQtWu8dUnKSi4jStK2HnoIvvUteO89uP56Z7IktZlhS5KafPWr8NvfRjNaf/kLDBkSd0WS2gHDlqSOLQyjP506RbNYhx4KV14Jnf3Xo6Tk8N8mkjquf/0LvvY1OP30aOnwi1+MuyJJ7ZAN8pI6nk2b4Lbb4JOfhGefhdzcuCuS1I45syWpY3n99Whz0r//HT77WbjrLujfP+6qJLVjhi1JHctbb0X7Zd19N1xyCQRB3BVJaucMW5LavyVL4Lnnov6s0aOjwLXXXnFXJamDsGdLUvu1YQP88IfRVg5lZbB2bTRu0JKURoYtSe3T3LlQVBSFrM9/HhYvhj33jLsqSR2Qy4iS2p///AdOPBH23hsefRTOOCPuiiR1YIYtSVmrorqW8soaltcluGjPF5j00p38s24dtw39BBu/th/LjxrAVz7ZibFxFyqpQzNsScpKFdW1TJq+mERDI19c/yTXP3QHnd9o4LAAbnl6LV/9/gBe5z+UzS4DYOzBRi5J8bBnS1JWKq+soaH+Q778wsNMuePXdF7awCYgJ4TcjSFFr6wDoL6xnqkLpsZbrKQOzZktSVlp+fvr+dMD13Ps20sID+0MQzrz4YwPyd0Y0tA5YN7AbluuXbluZYyVSuroDFuSssu//gV9+tBn7678ftjp3HX0WfzwiD/SL2cV1xy2PwNer2fewG68eGjXLb/Sq1uvGAuW1NG5jCgpOyQSMHkyFBbCHXdQWlzIrEEnUnXoMfyk8VzWh10Ys8+H3HfavtsFrbycPEqGlcRYuKSObpczW0EQ/A44HXgnDMNPbR7bB3gA6A+8CUwMw/D91JUpqcMKw2j7hu9+F958M9oz64wzGH9gXyDq3ZpRdzz75Hbh+5segP+8z9R992FlTkCvbr0pGVZic7ykWAVhGO78giA4AVgL3LtN2PoJ8F4YhrcEQXANsHcYhlfv6mZFRUXhvHnzklC2pA7jO9+BW2+FT34SfvUrOPnkuCuS1EEEQTA/DMOitn7PLme2wjD8exAE/T8yfCYwavP73wOzgF2GLUlqkTVrogOi99wTzj4bDj0UvvlNyM2NuzJJarXd7dnaPwzDFZvfrwT2T1I9kjqyMIQ//jHqyyori8ZGjYKSEoOWpKzV5gb5MFqH3OFaZBAElwVBMC8IgnnvvvtuW28nqb1auBA+8xm44ALo2xfOOSfuiiQpKXY3bP07CILeAJtf39nRhWEY3hmGYVEYhkU9e/bczdtJatduvx0+/WmoqYG77oLnn4djj427KklKit0NW48CF21+fxHwSHLKkdRhNDbC6tXR+5NOgm9/G159Fb7yFejkrjSS2o9d/hstCII/AXOAwiAIlgVB8GXgFuBzQRC8Bnx282dJapnZs+Hoo+HSS6PPhYUwdSrsvXe8dUlSCrTkacTzd/Cj0UmuRVJ7t3IlXH013Htv1JdVWhp3RZKUch7XIyk9/vpXmDABPvwQJk2Ca6+NtnaQpHbOsCUpNebMgVmzouXC0aNh6FAYOzY6cueww+KuTpLSxrAlKfnmzIl2eq+vjzYn/fvf4fjj4U9/irsySUo7H/mRlFzLl0fnGNbXbx2bNSu2ciQpbs5sSUqeOXOiJcMNG6Bz52hH+C5dojFJ6qAMW9Jumrl0JlMXTGXlupX06taLkmEljD14bNxlNW/RNKiaDKuXQY9+MPoGGDwxOd+9bl20P9bQodHGpF/7Glx+efTk4axZ0XE7I0Yk514fUVFdS3llDcvrEvQpyKe0uJDxQ/um5F6StLuC6LSd9CgqKgrnzZuXtvtJqTJz6UzKZpdR37h1qSwvJ4+y48oyL3AtmgYzLoeGxNax3Hw441dtC1wbNsCdd8J//Vc0i7V0aTSLlSYV1bVMmr6YREPjlrH83BxunjDIwCUpKYIgmB+GYVFbv8eeLWk3TF0wdbugBVDfWM/UBVNjqmgnqiZvH7Qg+lw1efe+r7ER/vAHGDgQvvOd6PXBB9MatADKK2u2C1oAiYZGyitr0lqHJO2Ky4jSbli5bmWrxmO1elnrxnflmWfgwgujZcMnn4QxY6InDtNseV2iVeOSFBdntqTd0Ktbr1aNx6pHv9aNN+dvf4O7747ejx4NTzwB8+ZBcXEsQQugT0F+q8YlKS6GLWk3lAwrIS8nb7uxvJw8SoaVxFTRToy+IerR2lZufjS+K9XVcOqpUZP7T34CGzdG4eqUU2I/LLq0uJD83JztxvJzcygtLoypIklqnmFL2g1jDx5L2XFl9O7Wm4CA3t16Z2ZzPERN8Gf8CnocAATR666a4998E847D4YNgxdegPJyWLgwaoTPEOOH9uXmCYPoW5BPAPQtyLc5XlJG8mlESVuFYTRztWhRtOP7d78LV14JPXrEXZkkpV2ynkbMnP+bKik+q1bBLbfA++/Db38LgwdDbS107x53ZZKU9VxGlDqytWujfbIOPhh+9rNoW4dNm6KfGbQkKSmc2ZI6qDl/uoWB37iRvVdv4LnhPdk4+WZOPOWbyfnyVO5Yvw13kJeUDQxbUkcxZw48/TQcdRQzj4D/t+pPXHdQF24f149Fh3Yl7z//TdnSg9re5P/RHetXvx19hqQGro/uIF9bl2DS9MUABi5JGcUGeakj+N//3XpAdKdOlEwextP96j92We9uvfnLOX9p271+8akoYH1UjwPge0va9t3bGHnL09Q2s4Fp34J8nrvm5KTdR1LH5XE9knatvh7uuAPOPDMKWpsdXN1MGCJJO+Ane8f6HXAHeUnZwrAltWe/+x18/euw//7R2YU5ObDHHiwdemCzlydlB/xk7FjfAu4gLylbGLak9mTFCrj6avif/4k+X3RR1Kf10kswaxZMmQJVVYw5/8bU7YDflh3rW8Ed5CVlCxvkpfbgtdeiXd5///voSJ3SUvjCF6BbNzjppOiaESOiP0BTC/zUBVNZuW4lvbr1omRYSXJ2wG9qgk/x04hNTfA+jSgp09kgL2W766+HH/0oWia85BK46io45JC4q5KkrGeDvNRRhSH89a/Rbu8ARUVwzTXReYa33WbQkqQMY9iSskVjIzz4IAwfDp/7HNx9dzQ+fnw0s9UrCc3tkqSkM2xJmS4M4a67YOBAmDgRPvgA7rwTvv3tuCuTJLWADfJSpvrwQ9hjDwgCePhh6NEjmtk666xoCwdJUlYwbHVgM5fOTM3TaNkoTWf5tcjKlTB1ajSbNW8e9O8P998fHQwdBPHUlKE8G1FSNjBsdVAzl86kbHYZ9Y3RkS0r1q2gbHYZQMcLXGk6y2+XXn8dfvpTuOceaGiAs8+O+rQA9torfXVkCc9GlJQt7NnqoKYumLolaDWpb6xn6oKpMVUUo6rJW4NWk4ZENJ5Kc+bAzTdHr3V1MGgQ/Pd/w8UXQ00NTJvmk4U7UV5ZsyVoNUk0NFJeWRNTRZLUPGe2OqgdnYGXlLPxsk2azvLbznPPbT0YOi8Pqqrg3nvh+OOhd+/U3bcd8WxESdnCma0Oakdn4CXlbLw2mLl0JmMeGsPg3w9mzENjmLl0Zupvmqaz/AB45x348Y/htFOiBvgwhPoE3P8r+PznDVqt4NmIkrKFYauDKhlWkrqz8XZTUx/ZinUrCAm39JGlPHCl6Sw/KiuhX79oA9KuCcgBAqLX+sqod0wtdtLAnq0al6S4uIzYQTU1wWfS04g76yNLaV2pOsvvvfeipcEDDoia3Y89Fr7zHch9FPLegbc3wpuN0D8HejdG94/rCcgs9Mwr77ZqXJLiYtjqwMYePDajnjyMtY9s8MTkBJ0whP/7P7j99qjBvb4eLr00Cls9esDPfgZlm3d+P6Bz9KdJKnvE2iF7tiRlC5cRlTEytY+sVS65BI47DqZPj94vXLj1WJ0m6ewRa8fs2ZKULQxbyhiZ2Ee2SwsWwDe+ES0ZApxzDtxxByxfDr/5DQwZ8vHfSVePWDtXWlxIfu72O+nn5+ZQWlwYU0WS1DyXEZUxMrGPrFnr1sEDD0RLhXPnQn5+dITOmDFw+um7/v1U9Yh1ME0bl7qDvKRMF4RhmLabFRUVhfPmzUvb/aSke++9aKPRujo48kj4+tfhggugoCDuyiRJSRYEwfwwDIva+j3ObEk7U18Pf/4zvPYalJXBPvvAVVfBCSdEG5B6VqEkaRcMW1JzXn0V7rwzOqdw1So44gi49lro0gWuuy7u6iRJWcSwJTWZMwdmzYI1a6IzCzt3hvHjo6XCk06CTj5PIklqPcOWOrbGRvjf/4Vf/hIefzz6nJsLl10WLRt6fI4kqY0MW+p4whBeeAHuvz/aeHT58ihgbdwY/Qygf3+DliQpKQxb6hjCMApVfTdvC3DhhfDmm3DaaXDeedCzZ7Rtw4YNUV/WqFFxVitJakcMW2rfXnkl2hPr/vth5crozx57RGMDBkRH6DSpqop6tkaNghEj4qpYktTOGLbUPv31r9EWDS++GG3PcOKJ8N3vwqZN0c+POurjvzNihCFLkpR0hi21D8uXw4MPwmc+A8OGQbdu0LUrTJ0aHaHTp0/cFUqSOijDlrLXu+9GG47efz/8/e9RX9aUKVHYGjECZs+Ou0JJkgxbyjIbN0b7X23aBIMGwb//DQMHwo03Ro3uhR5CLEnKLIYtZa6mTUaPOQbeeSeawXr1VXjppWiD0dtug4MPhsGDPTZHkpSxDFvKTHPmRLu2f/jh1rE+feDcc6PzCvPz4ayz4qtPkqQWMmwpM4Rh9OTgk09GIWrWLGhoiH4WBPCVr8Dtt3tkjiQp6xi2FJ/6enjkkShgPflktAcWwCc+Ee11tcceWzcZveQSg5YkKSsZtpQ+jY0wb14Usk48MWp2v/BC2HNPGDMGTjklem06JsdNRiVJ7YBhS6n1739DZSU88QT85S/w3ntw3HHw3HNRyKqujp4gzMn5+O+6yagkqR0wbCm5Ghqi3quioujzJZdEQWv//eGMM6LZq899buv1Rx4ZT52SJKWJYStmM5fOZOqCqaxct5Je3XpRMqyEsQePjbus1vnXv7bOXv31r7BmTTSjtd9+8MMfwk03wZAhreu5WjQNqibD6mXQox+MvgEGT0zd30FZqaK6lvLKGpbXJehTkE9pcSHjh/aNuyxJ2o5hK0Yzl86kbHYZ9Y31AKxYt4Ky2WUAmRm4mva9GjEi2vsqPx/uuw++9KXo5wccAOefH81ede8ejQ0f3vr7LJoGMy6HhkT0efXb0WcwcGmLiupaJk1fTKKhEYDaugSTpi8GMHBJyig+3hWjqQumbglaTeob65m6YGpMFe3Axo1wzz1RU/u110b7X910U/SzE06An/0s2mj0rbfgjjuirRvy83f/flWTtwatJg2JaFzarLyyZkvQapJoaKS8siamiiSpec5sxWjlupWtGk+blSujpcDDDoMPPoieDly/fuvPgwDq6qL3BxwAV1yR3PuvXta6cXVIy+sSrRqXpLg4sxWjXt16tWo8ZebNg6lToyXA/v2jcFVaGv1sr73g+9+Peq/y8qKnBvPy4ItfTF09Pfq1blwdUp+C5mdPdzQuSXExbMWoZFgJeTl5243l5eRRMqwkNTcMQ3j7bZg2LQpXTb7zHfjud6PtGI4+OloW/MEPtv78xhvhhhvg6adhypRo/6tUbskw+gbI/ch/MHPzo3Fps9LiQvJzt98yJD83h9JiDyOXlFmCMAzTdrOuA/LD439wCCUHn8XYUVPSdt9MlpanER9+GP74R/i//4Ply6OxvfeGd9+NZqoWLYJ994W+GdRU7NOIaoHrKxbzp+ffpjEMyQkCzj/mAP5r/KC4y5LUTgRBMD8Mw6I2f086w1b+gPzw0LJDydsUUjbAwJVUYQhLl0aBqunPk09GIepHP4Lf/Q6OPXbrnyFDIDc37qql3fbRpxEhmtm6ecIgn0aUlBRZHbYAejeG/OXSJWm7d2tl/P5XTzwBL7wQHW+zenV07M2770Y/69YtWg68446oyX3TprafK+hMkzLMyFuepraZZvi+Bfk8d83JMVQkqb1JVtiK7WnElRncLZZx+1+99x7893/DK6/Ayy/D4sXRU4KdOsGPfxzNWo0dG81YjRgBn/zk9sffJCNoue+VMkxzQWtn45IUl9jCVq9Ncd1513a2/1VKwtamTbBkSRSkmgLVK69Es1VXXAEbNsBVV0HPnjBwYHSW4Lx50e9t2AD//GcUxlJlZ/teGbYUk5wgoLGZmfmcIIihGknasVjCVt6mkJKDz4rj1i2Ssv2vVq3aPkwdcgh84xtRaBo+PApOQQADBkShqnfv6Pf23z9aIvzEJ6LPc+bA6NHR9V26wKhRbatrV9z3ShmouaC1s3FJiktaw1ZA1KuV6U8j9urWixXrVjQ7vkNNR9mccEL0VN8rr8CHH8KZZ0Y/LyqC+fO3Xr/tXlWdO0dPDPbrF/VYfXT39SDYGrQgWiqsqoruN2pUardhgKhHa/XbzY9LMelbkP//27v/2KrKO47j76+l0FIRYUTbASpMg+LGJqDZmE6dBpRGQWcUszlctumcm9UlGPcjiC5mS1yG6JxKNqMuTlQ20U0cuKnxBxatCv4cgtUwsNVNBEUqWvjuj+dcelvOvb3Qe+6vfl7JTc99zjm353x57umX53nOczKO2RIRKSUFHSA/efJkb2lpKdjv21s9x2xBmP/qqqN/wfS6SWGG9bY22Lo1TAT69NMhyers7P5Bhx4Ka9eG5QULYMeO0GJ1xBFw0EHdx1WVsp5jtiDMe3Xa9epGlKLR3YgikrSyHyBfctrbQ2LU1kZjWxvjXjuEda+v4PJvD6d+3wYWPjCIQ86b2X2f2lqYNSu0MO2ILvhmMHMmXHZZSKxSmhKaqLQQUgmV7kaUEpJKqK5dtoa3N3fwJ9/HgwAACmtJREFU2f1rmTNtnBItESk5lduy9cQTsGwZTJsGxx0Hq1eH6RLa2rpaptrbYcWK0EU3d26YHT2luhrq68MDlocMgfvvh1Wrwjiq+vrws6EhdBk2N3cfQ5X0DOsiIiKSuLKcZ2tyQ4O33HQTjB8PHR3hNXYsHHAAvPsuLF8OH3/cta6jA845J7QQrV4dHiOTvq6jA264ASZNgiVL4IILwv7btnW1NA0aBI8+GqZLuPDCkDilEqX6epg/PyyvWQPr13etGz48tFLlKjVmqxBjqERERCRx5dmN2N4OZ/S4C/G222D27NCFd955u+3y3AGfMOnwa2DzZnjyydB1V1MTftbVdSVEo0bBmWeG8hdegMcfD7Oqd3aGJKipKQxIr6uLP7Zx48Jrb7x4DzRfDds3QPOdUKcuNpFCWPLCRnUjikjJK/yYrX32gbPOCi1WtbUwYUIonziRx/55C795+Ua2VH3C9oH7sL3aGFi9nHmtU2g8vjE8jiaDB4e/w4KTW2n/qJ0TR9Tw25XVVH26o2tqhMGDkzmfYkz4qdncRVjywkbm3LuaT3eG1vmNmzuYc+9qACVcIlJSCtuNaOYttbUZxzRNXTw1dsqFhroGlp+1POPnxt09eHRrJ5dvO4bDz/xBst168z+fYVqE0XBZAo8j0p2BIgB86arlbO74dLfy/WurWXXl1CIckYhUmnx1I/bpOS5mdoqZrTGzdWZ2Ra87jBzJU3dczdSNVzLh9glMXTyVB1sf3LV6bycTjZvx/dmxA7hk3DPQ/P2QoCSl0BN+ZpvNXaQfiUu0spWLiBTLXidbZlYF3AicCowHzjWz8dn22TKslks/vpu2j9pwfNczB1MJV6ZJQ7NOJkqWJG1AVVe3XlIJV6aJPZOa8FOzuYuIiJSVvrRsHQOsc/dWd/8EWATMyLbDO9veyfjMQYCmiU3UVNV0W19TVUPTxOxzVGVM0jqjOxKTbPk5aW7oxktXXRvKk1Do5E6kRGW6WViPRhSRUtOXZGskkD5YaUNU1o2ZXWBmLWbW8unO+Ob9VMtU49hG5k2ZR0NdA4bRUNfAvCnzen34c2yStnMnTe9v7ipIquVnwtlhvNTQ0YCFn0mOnyp0cidSojINN9WjEUWk1CR+N6K7LwQWAuz3uf1iL4PpLVONYxt7Ta56Sm2/4PkFtG99m/rOHTS9v5nGj7Z1bZRky8+Esws3OF2zuYsAejaiiJSPviRbG4HRae9HRWUZHTj4QGqqanZ75mBv3YS52JWkZbpbr5JafgqZ3ImUqDnTxsU+G3HOtL2cL09EJCF96UZ8FjjMzMaY2UBgFvBAth2GDhq6V92Ee6TQ3XoiUhQzjxrJr878AiP3r8UILVp6CLWIlKI+zbNlZtOB64Aq4FZ3vybb9gV9NqKIiIhIH5TE43rcfSmwtK8HISIiIlKp+jSpqYiIiIhkp2RLREREJEFKtkREREQSVNAHUZvZh8Cagv3C8jEC+F+xD6LEKCbxFJd4iks8xWV3ikk8xSXeOHcf0tcPSXxS0x7W5GNUf6UxsxbFpTvFJJ7iEk9xiae47E4xiae4xDOzvEyhoG5EERERkQQp2RIRERFJUKGTrYUF/n3lQnHZnWIST3GJp7jEU1x2p5jEU1zi5SUuBR0gLyIiItLfqBtRREREJEF5SbbM7BQzW2Nm68zsipj1g8zs7mj9SjM7JG3dT6PyNWY2LR/HUypyiMtPzOxVM3vRzP5lZgenrdthZquiV9YHfJebHOJyvpn9N+38v5e2braZrY1eswt75MnKIS7z02LyupltTltXkfXFzG41s3fN7OUM683Mro9i9qKZTUxbV8l1pbe4fDOKx0tmtsLMvpi27q2ofFW+7rQqBTnE5AQz25L2PZmbti7rd6+c5RCXOWkxeTm6lgyP1lVkXQEws9Fm9mj0N/gVM2uK2SZ/1xd379OL8BDqN4CxwEBgNTC+xzY/BG6OlmcBd0fL46PtBwFjos+p6usxlcIrx7icCAyOli9KxSV6v7XY51DEuJwP/C5m3+FAa/RzWLQ8rNjnVKi49Nj+x4SHv1d6ffkaMBF4OcP66cBDgAFfBlZWel3JMS5TUucLnJqKS/T+LWBEsc+hCDE5Afh7TPkefffK7dVbXHpsexrwSKXXlejcGoCJ0fIQ4PWYv0V5u77ko2XrGGCdu7e6+yfAImBGj21mALdHy4uBk8zMovJF7r7d3d8E1kWfVwl6jYu7P+ru26K3zcCoAh9jMeRSXzKZBjzs7pvc/X3gYeCUhI6z0PY0LucCdxXkyIrI3R8HNmXZZAZwhwfNwP5m1kBl15Ve4+LuK6Lzhn5ybcmhrmTSl2tSydvDuPSL6wqAu7e5+/PR8ofAa8DIHpvl7fqSj2RrJPCftPcbYg541zbu3glsAT6T477lak/P7buEDDqlxsxazKzZzGYmcYBFkmtcvhE12y42s9F7uG85yvncou7mMcAjacWVWl96kylulVxX9lTPa4sDy83sOTO7oEjHVCxfMbPVZvaQmR0ZlamuAGY2mJAw/CWtuF/UFQtDm44CVvZYlbfrS6FnkJcYZvYtYDJwfFrxwe6+0czGAo+Y2Uvu/kZxjrDg/gbc5e7bzexCQqvo14t8TKVkFrDY3XeklfXn+iIZmNmJhGTr2LTiY6O6cgDwsJn9O2r9qHTPE74nW81sOrAEOKzIx1RKTgOecvf0VrCKrytmti8hwbzU3T9I6vfko2VrIzA67f2oqCx2GzMbAAwF3stx33KV07mZ2cnAz4HT3X17qtzdN0Y/W4HHCFl3Jeg1Lu7+Xlos/gBMynXfMrYn5zaLHk39FVxfepMpbpVcV3JiZhMI358Z7v5eqjytrrwL3EflDN3Iyt0/cPet0fJSoNrMRqC6kpLtulKRdcXMqgmJ1p3u/teYTfJ3fcnDILMBhMFhY+gaXHhkj20upvsA+Xui5SPpPkC+lcoZIJ9LXI4iDMw8rEf5MGBQtDwCWEuFDNjMMS4NactnAM3R8nDgzSg+w6Ll4cU+p0LFJdrucMKgVesP9SU6p0PIPOi5ke4DWJ+p9LqSY1wOIoyBndKjvA4Ykra8Ajil2OdSoJjUp743hKRhfVRvcvrulfMrW1yi9UMJ47rq+lFdMeAO4Los2+Tt+tLnbkR37zSzHwHLCHd13Orur5jZ1UCLuz8A/BH4k5mti/5BZ0X7vmJm9wCvAp3Axd69a6Rs5RiXa4F9gXvD/QKsd/fTgSOAW8xsJ6H18dfu/mpRTiTPcozLJWZ2OqFObCLcnYi7bzKzXwLPRh93tXdv8i5bOcYFwndnkUff+EjF1hczu4twF9kIM9sAXAlUA7j7zcBSwh1D64BtwHeidRVbVyCnuMwljIv9fXRt6fTwkOEDgfuisgHAn939HwU/gQTkEJOzgIvMrBPoAGZF36PY714RTiEROcQFwn9ql7v7R2m7VmxdiXwVOA94ycxWRWU/I/xHJe/XF80gLyIiIpIgzSAvIiIikiAlWyIiIiIJUrIlIiIikiAlWyIiIiIJUrIlIiIikiAlWyIiIiIJUrIlIiIikiAlWyIiIiIJ+j/h+2HRKDwMKAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,8))\n",
    "# Plot the graph\n",
    "ax.plot(X2,y2,'o',label='data')\n",
    "ax.plot(x2,y_fitted_new, 'r--.', label='OLS')\n",
    "# Set the legend\n",
    "ax.legend(loc='best')\n",
    "ax.axis((0,2,-1,50)) ## Set the area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dummy variable 哑变量\n",
    "$\\large \\beta_0x_0 + \\beta_1x_1 + \\beta_2x_2 + \\beta_3 x_3$\n",
    "\n",
    "Set the formula as dummy variable  \n",
    "$\\large Y = 10 + X + Z_1 + 3 Z_2 + 8 Z_3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample3 = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
      " 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n",
      " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]\n"
     ]
    }
   ],
   "source": [
    "groups = np.zeros(nsample3,int)\n",
    "groups[20:40] = 1\n",
    "groups[40:] = 2\n",
    "\n",
    "print(groups)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [],
   "source": [
    "dummy = sm.categorical(groups, drop = True)\n",
    "x3 = np.linspace(0,40, nsample3)\n",
    "X3 = np.column_stack((x3,dummy))\n",
    "X3 = sm.add_constant(X3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>1.0</td>\n",
       "      <td>38.383838</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>96</th>\n",
       "      <td>1.0</td>\n",
       "      <td>38.787879</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>97</th>\n",
       "      <td>1.0</td>\n",
       "      <td>39.191919</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>98</th>\n",
       "      <td>1.0</td>\n",
       "      <td>39.595960</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99</th>\n",
       "      <td>1.0</td>\n",
       "      <td>40.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      0          1    2    3    4\n",
       "95  1.0  38.383838  0.0  0.0  1.0\n",
       "96  1.0  38.787879  0.0  0.0  1.0\n",
       "97  1.0  39.191919  0.0  0.0  1.0\n",
       "98  1.0  39.595960  0.0  0.0  1.0\n",
       "99  1.0  40.000000  0.0  0.0  1.0"
      ]
     },
     "execution_count": 179,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(X3).tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [],
   "source": [
    "## beta\n",
    "beta3 = [10,1,1,3,8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [],
   "source": [
    "e3 = np.random.normal(size=nsample3) # error item by normal distribution\n",
    "y3 = np.dot(X3,beta3) + e3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>55.895872</td>\n",
       "      <td>1.0</td>\n",
       "      <td>38.383838</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-0.487966</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>96</th>\n",
       "      <td>57.087347</td>\n",
       "      <td>1.0</td>\n",
       "      <td>38.787879</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.299468</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>97</th>\n",
       "      <td>57.619158</td>\n",
       "      <td>1.0</td>\n",
       "      <td>39.191919</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.427239</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>98</th>\n",
       "      <td>59.241150</td>\n",
       "      <td>1.0</td>\n",
       "      <td>39.595960</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.645190</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99</th>\n",
       "      <td>59.502535</td>\n",
       "      <td>1.0</td>\n",
       "      <td>40.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.502535</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            0    1          2    3    4    5         6\n",
       "95  55.895872  1.0  38.383838  0.0  0.0  1.0 -0.487966\n",
       "96  57.087347  1.0  38.787879  0.0  0.0  1.0  0.299468\n",
       "97  57.619158  1.0  39.191919  0.0  0.0  1.0  0.427239\n",
       "98  59.241150  1.0  39.595960  0.0  0.0  1.0  1.645190\n",
       "99  59.502535  1.0  40.000000  0.0  0.0  1.0  1.502535"
      ]
     },
     "execution_count": 182,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_table = np.column_stack((y3,X3,e3))\n",
    "pd.DataFrame(result_table).tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[9.92971043 1.04239985 1.08376853 1.96273338 6.88320852]\n"
     ]
    }
   ],
   "source": [
    "result3 = sm.OLS(y3,X3).fit()\n",
    "print(result3.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.996\n",
      "Model:                            OLS   Adj. R-squared:                  0.996\n",
      "Method:                 Least Squares   F-statistic:                     8073.\n",
      "Date:                Sun, 26 Apr 2020   Prob (F-statistic):          3.32e-115\n",
      "Time:                        23:13:53   Log-Likelihood:                -132.82\n",
      "No. Observations:                 100   AIC:                             273.6\n",
      "Df Residuals:                      96   BIC:                             284.1\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          9.9297      0.199     50.015      0.000       9.536      10.324\n",
      "x1             1.0424      0.017     62.812      0.000       1.009       1.075\n",
      "x2             1.0838      0.205      5.287      0.000       0.677       1.491\n",
      "x3             1.9627      0.168     11.665      0.000       1.629       2.297\n",
      "x4             6.8832      0.307     22.414      0.000       6.274       7.493\n",
      "==============================================================================\n",
      "Omnibus:                        4.906   Durbin-Watson:                   2.066\n",
      "Prob(Omnibus):                  0.086   Jarque-Bera (JB):                4.561\n",
      "Skew:                           0.357   Prob(JB):                        0.102\n",
      "Kurtosis:                       3.766   Cond. No.                     4.08e+17\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The smallest eigenvalue is 3.23e-31. This might indicate that there are\n",
      "strong multicollinearity problems or that the design matrix is singular.\n"
     ]
    }
   ],
   "source": [
    "print(result3.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 20, -1, 40)"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAHWCAYAAABAA0zqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde5zVVb3/8deaCzCaNl7IYBTFG4KhjI4mUYmijZkHEc1LlpdMLMvqV2cUylOd4wVO2O10sfCSWpamIXnrkAdEU8eSmyISimbEoInmKOQAc1m/P74zwMAeZobZ39l7Zl7Px2Mes/d3vnvvNd8Q3q31+X5WiDEiSZKkdBTkegCSJEm9mWFLkiQpRYYtSZKkFBm2JEmSUmTYkiRJSpFhS5IkKUUdDlshhMIQwqIQwv3Nz4eGEP4UQlgRQrgzhNAvvWFKkiT1TJ2Z2foSsGyL5/8NfC/GeCDwJnBRNgcmSZLUG3QobIUQ9gY+BtzY/DwAxwN3N59yKzAhjQFKkiT1ZB2d2fo+cDnQ1Px8D6A2xtjQ/HwVUJblsUmSJPV4Re2dEEI4BXgtxrgghDC2sx8QQpgETALYeeedjzzkkEM6PUhJkqTutmDBgtdjjAO7+j7thi1gDDA+hHAyMADYFfgBUBpCKGqe3dobqMn04hjjDGAGQEVFRZw/f35XxyxJkpS6EMLfsvE+7S4jxhinxBj3jjHuB5wNzI0xngs8DJzRfNr5wO+yMSBJkqTepCt9tq4AvhJCWEFSw3VTdoYkSZLUe3RkGXGTGOM8YF7z45eAo7M/JEmSpN6jU2ErDfX19axatYr169fneiipGTBgAHvvvTfFxcW5HookSepmOQ9bq1atYpdddmG//fYjad/Vu8QYeeONN1i1ahVDhw7N9XAkSVI3y/neiOvXr2ePPfbolUELIITAHnvs0atn7iRJUttyHraAXhu0WvT230+SJLUtL8JWPvnWt77Fdddd1+bPZ82axXPPPdeNI5IkST1ZjwtbsxbVMGbaXIZOfoAx0+Yya1HGXqrpfb5hS5IkdUKPCluzFtUwZeYSamrriEBNbR1TZi7pcuC65pprOPjgg/ngBz/I8uXLAbjhhhs46qijOPzwwzn99NN55513eOKJJ7j33nupqqpi1KhRvPjiixnPkyRJatGjwtb02cupq29sdayuvpHps5fv8HsuWLCAO+64g8WLF/Pggw/y1FNPATBx4kSeeuopnn76aYYPH85NN93EBz7wAcaPH8/06dNZvHgxBxxwQMbzJEmSWuS89UNnrK6t69TxjvjjH//Iaaedxk477QTA+PHjAXj22We58sorqa2tZd26dVRWVmZ8fUfPkyRJfVOPmtkaXFrSqeNdccEFF/CjH/2IJUuW8M1vfrPN1g0dPU+SJPVNPSpsVVUOo6S4sNWxkuJCqiqH7fB7fvjDH2bWrFnU1dWxdu1a7rvvPgDWrl3LoEGDqK+v5/bbb990/i677MLatWs3PW/rPEmSJOhhYWtCeRlTJ46krLSEAJSVljB14kgmlJft8HseccQRnHXWWRx++OF89KMf5aijjgLgqquu4v3vfz9jxozhkEMO2XT+2WefzfTp0ykvL+fFF19s8zxJkiSAEGPstg+rqKiI8+fPb3Vs2bJlDB8+vNvGkCt95feUJKm3CCEsiDFWdPV9etTMliRJUk9j2JIkSUqRYUuSJClFhi1JkqQUGbYkSZJSZNiSJElKkWELWLVqFaeeeioHHXQQBxxwAF/60pfYuHEj8+bN45RTTtnm/Pvvv5/y8nIOP/xwRowYwc9+9rMcjFqSJGXdm29m/S37fNiKMTJx4kQmTJjACy+8wPPPP8+6dev4+te/nvH8+vp6Jk2axH333cfTTz/NokWLGDt2bPcOWpIkZUd1NVx9NUybBiefDIMHw5o1Wf2IHrUR9SbV1TBvHowdC6NHd+mt5s6dy4ABA7jwwgsBKCws5Hvf+x5Dhw7luOOO2+b8tWvX0tDQwB577AFA//79GTZsx7cLkiRJOXL//XDaadDQkDzfc0+oqoKC7M5F5V/YyjRLdOaZcOml8M47MGYMPPMMNDUlF+Oww+BLX4ILLoDXX4czzmj92nnztvtxS5cu5cgjj2x1bNddd2XIkCGsWLFim/N33313xo8fz7777su4ceM45ZRTOOeccyjI8v8wkiQpBRs2wD//CYMGJRmhJWgVFMCXvwxtrGx1Rc9LCG+9lQQtSL6/9Va3D+HGG29kzpw5HH300Vx33XV8+tOf7vYxSJKUtlmLahgzbS5DJz/AmGlzmbWoJtdD2nHPP5/MWu29N3z2s8mx00+HAQOgsBD694fjj0/lo/NvZmt7M1E77QS33w7jxsHGjdCvX/K8ZSlxzz3bncna2ogRI7j77rtbHXv77bdZuXIlBx54IH/4wx8yvm7kyJGMHDmST33qUwwdOpRbbrmlU58rSVI+m7Wohikzl1BX3whATW0dU2YuAWBCeVkuh9a+LcuNXn8dvvMdeOQRKCqC8ePhkkuS80aPhrlzs1aa1Jb8C1vtGT0a5szJ2oUZN24ckydP5rbbbuO8886jsbGRr371q1xwwQXstNNO25y/bt065s+fv6kofvHixey7775dGoMkSflm+uzlm4JWi7r6RqbPXt7lsDVrUQ3TZy9ndW0dg0tLqKoclr0AV12dzFDV1yeTMueeC3//O0ydmpQcvfe9rc8fPTq1kNWi54UtyOqFCSFwzz33cOmll3LVVVfR1NTEySefzLXXXkt1dTVz5sxh77333nT+r3/9a7797W9zySWXUFJSws477+ysliSp11ldW9ep4x2V2ozZO+/AXXfBN74B69cnxzZuhCFD4IUXMha9pxr6ttAzw1aW7bPPPtx3333bHB87dix1ddv+ofrQhz7UHcOSJClnBpeWUJMhWA0uLenS+2Z9xuxf/4IrroBf/jKp495nHyguTuq6+/WDE05oM2h11zJpzyuQlyRJqauqHEZJcWGrYyXFhVRVZm531NFi+qzMmK1blywXQlLP/cc/wimnJHVZf/tb8v2qq5KyozZWwrYX+rLNmS1JkrSNltmdjiyzdWaWaIdnzKqrk5viVq1KQlRREaxeDSUlsHBhckdhiw6UG6W1TJqJYUuSJGU0obysQ0tqnVkarKoc1iqYwfZnzAD48Y/hsssgxuT5ySfDlVcmbRugddDqoLSWSTPJi2XE2HLxeqne/vtJkvq2zswSTSgvY+rEkZSVlhCAstISpk4c2TqUxQhPPpkUtgMsXbo5aBUWwgc/mMxchbDDY+7sMmlX5Hxma8CAAbzxxhvssccehC5ctHwVY+SNN95gQEv6liSpl+nsLFGbM2Zvvgm/+AXccAM8+yx87nPwk5/AJz8Jt9yyucfmdvYk7ugdhp1ZJu2q0J2zLhUVFXH+/PmtjtXX17Nq1SrWt9ym2QsNGDCAvffem+Li4lwPRZKkrNu6ZguSWaJtZqy2tPU+x5ddBjfemLRtqKiASZPg7LNhl10yn5+tcWxHCGFBjLGi0y/c+n1yHbYkSVLP16meVS2NRzdsSOqu5syB++6D2lq4+GIoL9+hMYyZNjfjDFtZaQmPT+78VjzZCls5X0aUJEk9X4eK6ZuaktmpL3+5dePRefPg2mvbfFlHg1x33mHYGYYtSZKUvuXLk15YK1YkS4NFRUnRewdqsFJvK5GyvLgbUZIk9TJNTTB7NsycmTzfbz8YPjwpgP/HP+DRR9ttPAqdaz7anXcYdoYzW5IkKTuqq+Hee+GNN+Chh+Dll+Goo2DiROjfP/lZiw7uc9zZthLQPXcYdoZhS5KkHqy7NlNuV3U1HHss1NcnzysqYNo0mDChS2+btbYSOeQyoiRJPVRLPVNNbR2RzfVMbe1LmHUrV8K3vpXMYM2bB43Ny32Fhcls1llnJTNaXZCvS4Od4cyWJEk9VGe2ycma+np44IGk8ejvf58c22+/pMi9f/8ONR7tjHxdGuwMw5YkST1Ut7Q62LKZ6JFHwsEHw9/+BoMHJ/sTfvrTSdiCpNi9ncajOyIflwY7w7AlSVIPlXqrg0cfhRNPTGazWpqPXnppclfhRz+atG/YUgeL3vuadmu2QggDQgh/DiE8HUJYGkL4z+bjt4QQ/hpCWNz8NSr94UqSpBap1TM9/zxcfjmcfHKyLBjj5uajl18O//Zv2wYttakjV2oDcHyMcV0IoRh4LITQvEhLVYzx7vSGJ0mS2pJKPdPdd8PHP54UuY8ZA3/6EzQ0ZLUOq69pN2zFZPPEdc1Pi5u/um9DRUmS1KYu1zMtW5YUux91FJxzDowbl2ydc8EFMGhQhzaA1vZ1qPVDCKEwhLAYeA14KMb4p+YfXRNCeCaE8L0QQtfu7ZQkSd3j4YeT2avDDoMRI+BHP0pCF8Buu8GUKUnQgiRgTZli0OqCDi24xhgbgVEhhFLgnhDC+4ApwKtAP2AGcAXwX1u/NoQwCZgEMGTIkCwNW5Ik7ZDqajjhhGQ7nRDg85+Hb3wD3vOeXI+s1+pUU9MYYy3wMHBSjPGVmNgA/Bw4uo3XzIgxVsQYKwYOHNj1EUuSpI5btw5uuinp7v7mm8mSYIuCAigrM2ilrCN3Iw5sntEihFACnAj8JYQwqPlYACYAz6Y5UEmS1AkLFsBnP5v0w/rMZ+D115P+WC3NRwsLLXrvJh1ZRhwE3BpCKCQJZ7+JMd4fQpgbQhgIBGAx8NkUxylJktqydRH7Sy8lexOWlMCZZ8LFF8MHPpAsG0JqzUeVWUfuRnwGKM9w/PhURiRJkjruiSfg+OOTPlgFBfDHPyYB6q67ktqs0tJtX2Pz0W7lRtSSJPVEb76Z3EU4cSJs2JA0Hm1qSu40BDjjjMxBS93OsCVJUk/REqgAvvtduOyyJFAVFyc1WAMGwHHH5XaM2oZhS5KkfPfGG/C978Ghh8L99yfHLr00KYL/y1/gkUfgqquSWqx2lgdnLaphzLS5DJ38AGOmzWXWoppu+AX6Njc2kiQpHz3xBPz858kdhI88ktRkHXMM7LRT8vNBg1o3Hu1ADdasRTVMmbmEuvpGAGpq65gycwlA17rQa7sMW5IkdYNZi2o6tofhhg2wcGFS3F5Xlxw744yk8ejIkV0aw/TZyzcFrRZ19Y1Mn73csJUilxElSUpZy4xSTW0dkc0zSpuW8Jqa4A9/SLbQ2X9/nvvFTBrWbwCgIRSw9L0HdjloAayurevUcWWHYUuSpJS1NaN0493VcM01cMABUFkJDz/MiuNP4TtvlVJfWERDKKC+sIir170nK7VVg0tLOnVc2WHYkiQpZS0zR0fULOPzT9zJmL8uBmCnl1+EK6+E/feHX/8aamo4/9CzmLPPKM49+xq++6FPcu7Z11C918FMn728y+OoqhxGSXFhq2MlxYVUVQ7r8nurbdZsSZKUssGlJey1dCF3/GoKxU0NNIYCzjz3v6l5XwW8+GIStpq1BLOFZcNZWDZ8m+Nd0VKX1aHaMWWNYUuSpJRVVQ7j5dm3UtzUQABCjHywZin7f+Uc2L910BlcWkJNhmCVraW+CeVlhqtu5jKiJEkpm1BexpgTKghAI4H64n4ced5pGUOPS329jzNbkiR1g6P23hWAwi9/icIzz+TYNvpiudTX+xi2JEnqDg0NMGIEXHddsrXOdrjU17u4jChJUne46CJYurTdoKXex7AlSZKUIsOWJElpe+452HdfmDcv1yNRDhi2JElK2/z5sHIlvOc9uR6JcsCwJUlS2hYsgJ13hmG2b+iLDFuSJKVtwQIYNcri+D7KsCVJUpoaG2HRIjjyyFyPRDli2JIkKU3vvAMXXAAnnZTrkShHbGoqSVKadtkFfvzjXI9COeTMliRJaXrttaR7vPosZ7YkSUrTGWckhfEPP5zrkShHnNmSJCktTU1Jcfz73pfrkSiHDFuSJKXl+edh3TrvROzjDFuSJKVl/vzku2GrTzNsSZKUlgULoKQEhg/P9UiUQxbIS5KUlrPOSuq1ivznti/zf31JktJyzDHJl/o0lxElSUrDmjXwf/+XdJBXn2bYkiQpDQ89BCeeCCtW5HokyjHDliRJaViwAAYMgBEjcj0S5ZhhS5KkNCxYAIcdZnG8DFuSJGVdUxMsXGh/LQGGLUmSsm/FCli71rAlwNYPkiRl39ChyTLiPvvkeiTKA4YtSZKyrbgYjjgi16NQnnAZUZKkbPvxj2H27FyPQnnCsCVJUjY1NcHXvgazZuV6JMoThi1JkrLpxRfh7bctjtcmhi1JkrJpwYLku2FLzQxbkiRl04IF0K8fHHporkeiPGHYkiQpm156Kekc369frkeiPNFu64cQwgDgUaB/8/l3xxi/GUIYCtwB7AEsAD4VY9yY5mAlScp7v/1t0tBUataRma0NwPExxsOBUcBJIYRjgP8GvhdjPBB4E7govWFKktSD7LJLrkegPNJu2IqJdc1Pi5u/InA8cHfz8VuBCamMUJKknmL2bDjnHFizJtcjUR7pUM1WCKEwhLAYeA14CHgRqI0xNjSfsgooS2eIkiT1EHPnwsyZsOuuuR6J8kiHwlaMsTHGOArYGzgaOKSjHxBCmBRCmB9CmL/GpC9J6s0WLICRI6F//1yPRHmkU3cjxhhrgYeB0UBpCKGlwH5voKaN18yIMVbEGCsGDhzYpcFKkpS3YoSFC+2vpW20G7ZCCANDCKXNj0uAE4FlJKHrjObTzgd+l9YgJUnKe3/9K7z5pmFL22i39QMwCLg1hFBIEs5+E2O8P4TwHHBHCOFqYBFwU4rjlCQpv/3zn3D44XDUUbkeifJMu2ErxvgMUJ7h+Esk9VuSJKmiAhYvzvUolIfsIC9JUjbEmOsRKE8ZtiRJ6qoY4cAD4brrcj0S5SHDliRJXfXyy8meiO96V65Hojxk2JIkqasWLky+eyeiMjBsSZLUVQsWQFFR0tBU2ophS5KkrlqwAN73PhgwINcjUR7qSJ8tSZK0PePGuUWP2mTYkiT1arMW1TB99nJW19YxuLSEqsphTCgvy+6HXH55dt9PvYphS5LUa81aVMOUmUuoq28EoKa2jikzlwBkL3D985/Qr593IqpN1mxJknqt6bOXbwpaLerqG5k+e3n2PuQ734E994SNG7P3nupVDFuSpF5rdW1dp47vkAUL4JBDktktKQPDliSp1xpcWtKp450WIzz5ZNL2obo6O++pXsewJUnqtaoqh1FSXNjqWElxIVWVw7r+5u+8A2efDW+9lTQ1HTfOwKWMDFuSpF5rQnkZUyeOpKy0hACUlZYwdeLIHS+Or6uD+fOTxyUlMG9e8jjGpGar5bm0Be9GlCT1ahPKy7p+5+GSJXDDDfCLX0AIUFOThK2774bKyiRo9esHY8dmZczqXQxbkiRtqbo6maEaOxYaGpIeWk8+mYSp00+HSZM2d4r/0IdgzpzN548enbtxK28ZtiRJalFdDccdB/X1SUf4H/wgqcn67nfhU59KWjxsbfRoQ5a2y7AlSdLatfDrX8N//ids2JAc27gR1qyBpUuTpUNpB1kgL0nq2/7932HQILjkkmR5sLgYCguTZcPjjjNoqcuc2ZIk9Thd2u+wthZ+9zs477wkSBUWJi0cLr4Yjj46qc+yBktZFGKM3fZhFRUVcX7LLbOSJO2Arfc7hKR31nZbOjzxBNx2G6xcmQSpujp46imoqOieQatHCiEsiDF2+Q+JM1uSpB5le/sdZgxbv/lNMnPVMrlw6qnwH/8BRx7Z5bF0aYZNfYZhS5LUo7S732GM8OijyV2E48fD8uWbg1ZhIbz//VkLWlvOsNXU1jFl5hIAA5dasUBektSjtLWv4Yii9TB9erIp9Nix8K1vJT844YSkAWlL0XuWGo9ub4ZN2pIzW5KkHqWqchhTZi5h+MvPcszKJTw5ZCSjX/kLX33kNmiohw9+EL7+dTjjjOQFo0en0ni03Rk2qZlhS5LUo0woL2OvRx6i4trJhNhEQ1E/nr9sMgWHfQE+8xkYMWLbF6XQeHRwaQk1GYJVWzNv6rtcRpQk9QyNjfD738NppzH6K5+muKmRohgZ0NTAYXv2T7q8ZwpaKamqHEZJcWGrYyXFhVRVDuu2MahncGZLkvq4HnFH3caNcOihsGIFDBwIn/gE/Pa3ybY6OdoAuuUa5f21U84ZtiSpD8vbO+oaGuDBB5MGo9demwSqCy+Egw9O7jDs1w8+//mcNx+dUF5muFK7bGoqSX3YmGlzM9YdlZWW8Pjk47t3MNXVcM89yX6Ef/gDrF6dbKPz3HNQWtq9Y5GwqakkKQvy5o666upkhmrjxuT56NHwk5/Axz4GRf5TpZ7NP8GS1Ifl9I66FSvgxhuhvJyljy3mkPoGCoGGUMDyI4/l0FNPTX8MUjfwbkRJ6sO6/Y66DRvgzjuTRqMHHQTXXcfzs//I1evew8bCIhpCAfWFRVy97j3MWlSTzhikbubMliT1YZ29o67Tdy5WV7cuYh8/PqnH2ndfuPpquPBCLrztL9TU1nHu2ddsalK6cK+DWdnWXodSD2OBvCSpQ7a+cxGSWbCpE0dmDkXz5kFlZVKHVVKSdHH/17+gqSmZ2SpIFleGTn6ATP8SBeCv0z6Wyu8idUS2CuRdRpQkdUiH9wJ89ln44hdpqDxpU8F74/oNLP3VvUnI+shHNgUtaLs+zE7s6i0MW5KkDunQnYsvvwwjR9L405/xp8HD2VBYTEMoYON26rDyqRP7rEU1jJk2l6GTH2DMtLnWjSkrrNmSJHVIpjsXD/3Hi1y0bA5cPAtuuAH22w9uv51/W1bCc/X9OKJmWbt1WPnSiT1vG7yqx7NmS5LUIS1hpPz5BZy36H4OeGMVB/1zFY39+1N47rlJG4cQgJ5Zh5VXDV6VF2xqKknKKJW9DmNkwqjB7Pb0AsZc+02KmhppIrDi3M9w4A+/Dbvt1ur0nPbv2kF50+BVvY41W5LUi7TMPtXU1hHZvBS2w7VHtbXw4x9DeTncfz/HvrKUomTyioLCAg48dP9tghbkVx1WR1mor7QYtiSpF+nwHYPbEyM88USy8fPgwfCFLyR3DxYXJ/2y+vWDwsLk+9ixGd9iQnkZUyeOpKy0hECyFNdmi4g80RMDonoGlxElaTtSWZJL0Q4vhVVXJ32wxo2D978fzj0XXn8dPvUpmDQJjjxy87lz5rRuVNqGCeVleX2ttpYvhfrqfSyQl6Q2dLqJZx7odJF3jHD99fDFL0Jj4+bmoyUlcOCB8K53dcOopfzUbU1NQwj7hBAeDiE8F0JYGkL4UvPxb4UQakIIi5u/Tu7qYCQpn2RlSa6bdXgp7I034LrrYPhw+Pznk6AFSRPSefNg1CiDlpQlHanZagC+GmMcARwDfD6EMKL5Z9+LMY5q/nowtVFKUg70xLvTtlsr1dQE69YlJy5dClVVsOeecOWVyUxWO3VYknZMuzVbMcZXgFeaH68NISwD8nP+XJKyqCe2L4AMtVKvvgpTpyZ9sE4+GX74Q/jQh2DZMjjkkOSck0/uUB2WpM7r1N2IIYT9gHLgT82HvhBCeCaEcHMIYdt7fyWpB+vRd6dVV8NFFyXhaZ994GtfS74fd1zy8xA2By1IAtaUKQYtKQUdvhsxhPAu4LfAl2OMb4cQrgeuAmLz9+8An87wuknAJIAhQ4ZkY8yS1C165N1pr7wCf/1rsuHz+vVJAfwnPgHf+AYM6wEhUeqFOhS2QgjFJEHr9hjjTIAY4z+2+PkNwP2ZXhtjnAHMgORuxK4OWJK6U49oX9DQAL//PcyYAQ8+CJdckhS6x5jUYb3vfQYtKYc6cjdiAG4ClsUYv7vF8UFbnHYa8Gz2hydJatNbbyUzVvvuC+PHw/z5cMUVcNJJHWo8Kql7dGRmawzwKWBJCGFx87GvAeeEEEaRLCO+DFySygglSYmWxqPDhsHHP54EqZ/8BI4+Gn70IzjllKTLO3S48aik9HXkbsTHSDZq35qtHiSpu9x1V1J71dCQFLeXlcEHPpDUZ+2yy7bnjx5tyJLyhHsjSlI+e+SRpNj9zDOToAVJ2Jo3L3mcKWhJyiuGLUnKN8uXw5o1yeM1a2DFCrj4YhgwIKnD6t9/cwsHSXnPsCVJ+WD9erj9djj22KT/1fXXJ8dPOw1efDG503DuXLjqqqQeyyVCqcfocJ8tSVJ2zVpUw4MzZnLxfT9lxJqX2XnjO3DAATBtGlxwQXJS4RZNVa3Dknokw5Ykdbd33uGJm3/LnQte5+ZfTqZ/w0ZiKOCqj3yWkdd+jQlH7rPNS2YtqulZzVUlbeIyoiR1l8WL4fOfh0GD+MBl53Hsc49R3NhAAUkPnf7r32H6Qy9s87JZi2qYMnMJNbV1RKCmto4pM5cwa1FNd/8GknaAM1uSlLZnnoHPfAaeeiopbv/4xzmr6X3Uh0LOX3g/NDZQX1jEk0NGsjrDxtfTZy+nrr6x1bG6+kamz17u7JbUAxi2JCnbnngiKXY//HCYNAkGDYLGRvjBD+CTn4Tdd2fVtLnU1NZx7tnXcMzKJTw5ZCQLy4ZTVlqyzdtlCmDbOy4pvxi2JClb3noLrr4avvOdZF/CggIYOTIpal+woNWpVZXDmDJzCQvLhrOwbDgAJcWFVFVuu4fh4NISajIEq8EZgpmk/GPNlqS8NGtRDWOmzWXo5AcYM21u/tcnXX01DB4M112XBC1o3Xx0KxPKy5g6cSRlpSUEoKy0hKkTR2ZcFqyqHEZJcWGrY20FM0n5x5ktSXmnpSC8pU6ppSAcyJ8apX/+M1kqPP982HXXZPucc8+FY46BL3wBNm5sdxPoCeVlHfp9Ws7xbkSpZwqx5f+BdYOKioo4f/78bvs8ST3TmOZ6pq2VlZbw+OTju39A1dXJDNWxxyZb5syYAXffDRs2wJ13JlvpZDrfTaClHi2EsCDGWNHV93FmS1LeyauC8OpqGDcumalqakqWCHfdFS66KNlCZ9SobV9j81FJWzBsSco7eVEQ3tQEDz+c3EG4cWNyN2EIcPrpcOutsPPOGV9m81FJW7NAXlLeyWlB+KuvJtvlHHwwnHBC0sahX79k25wBA+CrX91u0LL5qKStGbYk5Z3O3KmXVbfdBvvsA1OmwN57wy9/CatWJRs/d2AD6O01H5XUd7mMKI5KW30AACAASURBVCkvdfROvR1WXQ333gtvvJE0Gv3wh5Mg9eUvJ93eh20xi9bBGqy8qjWTlDcMW5L6loYG+P734YorkrosSIreP/xhOOggmD59h986L2rNJOUdlxEl9S0f+hBUVW0OWoWFsP/+WXlrm49KysSwJan3qq+He+5Jmo3W1yfHLrsMpk6FkpIkaLXTeLQzclZrJimv2dRUUu/R0kz0oINg4UL4+c+TuwvLypLi9i3rsGw8KqkdNjWVpC21NB/dsCFZIgwBTjkFJk2Ck06Coq3+urPxqKRuYtiS1LM9/zzccAMsXry5y3sIcPnlSb8sScoxa7Yk9Tzr18OvfpUsAQ4bltxdWFLSuvnoqafmepSSBDizJakn+vKX4Wc/S+4inDoVLrgA3vte67Ak5SUL5CXlr+pq+MMfkt5YDz8MP/whlJfD0qVJ4ftxx0GBE/SS0mGBvKTe7bbb4KKLkqAFyTY6b7yRPD700ORLknoAw5akbjNrUQ3TZy9ndW0dg0tLqKoc1roHVYxJcXtdHVx88eagVVAAn/1ssjF0d4xDkrLIsCVpG2mEkVmLapgyc8mmjZprauuYMnMJABOaXoUZM2DJEnj88aTY/brrki11Nm5MCt+PO67Lv1e74zBwSUqBYUtSK2mFkemzl1NX38gRNcs4ZuUSnnnvQexX+wrDb74MVq9IAtZZZ8E778DOOyed3isqsl7w3jKOLdXVNzJ99nLDlqRUGLYktZJWGFldW8cRq5Zx+51fp7ixgaZQQL+mBp57z1D48Y/hE5+A0tLWL0qh8ejqDBtFb++4JHWVYUtSK6mEkTff5IvP/S/nzv0l/Rs2UgA0ALeVn8zPPv4VHr903I6/d7OOLn0OLi2hJsPvMri0pMtjkKRMvGdaUitthY4dCiOPPQbnnQeDB/P/7vsRawe8i4aCIhpCAfWFRfz+8BOoOumQLo5489JnTW0dkc1Ln7MW1WxzblXlMEqKC1sdKykupKpy2DbnSlI2OLMlqZWqymGtaragk2Fk7VrYZZfk8TXXwBNPwIUXwsUXs4T38N8zZnLgc/NZMaKCsyZNzEqdVGeWPlueezeipO5iU1NJ2+j03YiPPw4//zm8/HISrv7yFxgyBP72N9hzz6TgPUVDJz9Apr/JAvDXaR9L9bMl9V42NZWUmgnlZR2b6amtha9/Ha6/PumRBfDxjyf7EwLsu296g9yCdViS8pk1W5I6p6kJXnsteVxXR9NPf0rLDHlDKGDpXgdAWfcuyVmHJSmfGbYkdczq1UkN1gEHwNlnAzDr1SY+f8Z/sL6o36ai96vXvSdjYXqaJpSXMXXiSMpKSwhAWWkJUyeOtA5LUl5wGVFSZtXVSUPRd7872Qz6/vuhsRHGjYNJk4CkyLxm6FH84+xrOGblEp4cMpKFex3Myhw0CO3w0qckdTPDlqRt3XMPnHtuslVOQQG8611QVZVsDH3ggZtOa+m9tbBsOAvLhm9zXJJk2JLUor4eHngAbrgBHnww2RC6pej9K1+BK6/c5iUWpktS+6zZkvq6+vrkjsIhQ+C002DxYrjgAujfP7mrsF+/ZOkwAwvTJal9zmxJfdHGjfDMM8lGz8XF8H//lzyeNAk++lEoKkoet7MJtA1CJal9NjWV+orqarj7bnj1VXjoIVi3Dl55JSmA37gxmcGSJG3SbU1NQwj7ALcBewERmBFj/EEIYXfgTmA/4GXgzBjjm10dkKQUVFfDsccmS4YAH/4wTJ6cFL6DQUuSUtSRmq0G4KsxxhHAMcDnQwgjgMnAnBjjQcCc5ueS8tEDD2wOWoWFcNJJyXJhYeH2XydJ6rJ2w1aM8ZUY48Lmx2uBZUAZcCpwa/NptwIT0hqkpC762MdaF7yPHZvrEUlSn9GpAvkQwn5AOfAnYK8Y4yvNP3qVZJlRUj7517/gscegshIefrjdgndJUvZ1uPVDCOFdwG+BL8cY397yZzGpss9YaR9CmBRCmB9CmL9mzZouDVZSJ8QIF14IH/sYf7jvCcY8UsfQtw5jzCN13b6djiT1ZR0KWyGEYpKgdXuMcWbz4X+EEAY1/3wQ8Fqm18YYZ8QYK2KMFQMHDszGmCV1xDXXwF138ewXJvOlP79NTW0dEaiprWPKzCUGLknqJu2GrRBCAG4ClsUYv7vFj+4Fzm9+fD7wu+wPT9IO+d3v4D/+Az75SS7Z6zjq6htb/biuvpHps5fnaHCS1Ld0ZGZrDPAp4PgQwuLmr5OBacCJIYQXgBOan0vKtdWr4ZOfhKOOghkzWP3W+synuX+hJHWLdgvkY4yPAaGNH2few0NS7gweDN//ftLeoaTE/QslKcfcG1HqLerrYXnz0uBFF0FZsmWO+xdKUm65N6LUQ81aVNNqT8JfLP4F+99/FyxbBvvss+k89y+UpNwybEk90KxFNUyZuWRT4fuH5t3D/rNv4YVPXcJBWwStFhPKywxXkpQjhi2pB5o+ezl19Y0cUbOM05fM4cxnHuKRoUdw5SGn88dcD06S1IphS+qBVtfWcUTNMn51x9fp37CRCNxw1ARWvb0x10OTJG3FsCXlka3rsDLWVjU1cfy6lQxbuYSixgYC0BgKOOzVFfy1fExOxi1Japt3I0p5oqUOq81O76+8AtdeCwceyE0/vpS/D9yH+sIiGkIB9YVFLNp/lHcYSlIecmZLyhMtdVhbqqtv5Fe/epgJ//VbuO8+aGyE446Da67hxKFH86WBe3Hgc/NZMaKCsyZNtAhekvKQYUvKEy0d3Y+oWcaJLzzJC3sMYebIcby4oQDmz4evfhU+8xk46CAAxgPjj7kshyOWJHWEYUvKE/vs0o+PPvQrLn/0VgpipCkEXt59MP849Aj429+gwFV/SeqJDFtSPvjZz5j9g29SsuYfRJL9sSKBD9YsZf+vnGPQkqQezL/BpU6ataiGMdPmMnTyA4yZNndzAXtnbNwId98Na9cmz5uaKDm6gme/OIUNxf1pCAU0FBVz5HmnWYclST1ciDF224dVVFTE+fPnd9vnSdm2ded2SPYZnDpx5PZDUXU1zJsHBx6Y1F/dcgu89hr8/OdwwQWZzx07FkaPzv4vIUnqkBDCghhjRVffx2VEqRPaumNw+uzlbYet6moYNw7qkgJ4Cgpg/Hi4+GKorNz2/NGjDVmS1IsYtqROaLljsEPHly2DhQtZ+uSzDFu/gSKgicCK8z/HwTf/KN2BSpLyhjVbUicMLi3Z/vG6OrjtNvjQh2DECOov+Szfrn33puajG4qK+WY8YMfqvCRJPZIzW1InVFUOy1izVVU5DH73u6T+qrY2qc369reZ8Nb+LG0YwLlnX8MxK5fw5JCRLNzrYFZub9lRktSrGLbUa3Von8FOann9gzNmMnzJk+xeDPtdcBbHlpfBzv+Cj34UJk2CY4+FEHhu8gMALCwbzsKy4Zvep63lSElS72PYUq+09V2DLfsMAl0PXE8/xIQb/l+ydQ7Anhvh/PFw8MHwq1+1OndwaQk1GYJVW8uRkqTex5ot9Urbu2uwS049FS68cHPQKiiA8vI2T6+qHEZJcWGrY5uWHSVJfYJhS71Sp+4abEuM8Oc/w//7f9DQkBz72MfgK1+BkhIoLIT+/ZONodswobyMqRNHUlZaQgDKSkva78klSepVXEZUr7RDy3ctzUSPPBJeeAFmzIBnnoGddoLzzktmsCZNSs4944wONx6dUF5muJKkPsywpV5pu3cNZvDILffy/klnUtywkYIYCQBHHAHXXw+f+ATsumvrF9h4VJLUQYYt9UotM0nt3o34xhssufaHxHsfpKihnsIYaQJurziFd834qTNSkqQuM2yp12pz+S5GeOQRuOEGuPtuRm7cyFNlw6kvLILGBuoLi5h5yLH8w15YkqQsMGyp77n8crjuOnj3u2HSJD76r2Ese89QjqhZtrnxaNlwgr2wJElZYNhS7/b443DzzfDyy/Cd78CoUXDuuTByZFLkvtNOvD1tLtTWbdN41F5YkqRssPWDeqdXXoHPfS7Zo/Dmm2HuXHgg6ebOqFHJ3YU77QTYC0uSlC5nttT7bNwIhx4Kb765+VhhYdKANIMOF9NLkrQDDFvq+f7+92T26k9/Smav+vVLit8bGpJu7xs3JsfGjm3zLeyFJUlKi2FLPUtL49EPfjCZuZoxA37/e2hqghNPhLfegtJSOP305PwhQzrcfFSSpDQYttRzVFfDuHHJTFVhYfJ90CCYPBkuugj233/b19h8VJKUY4Yt5b/6erj3XpgyBdavT/pkAXzyk/Dzn0ORf4wlSfnLf6XUIbMW1XR/AfmKFXDjjXDLLfCPf8DAgUmwampKarAuvdSgJUnKe/5LpXbNWlTTap/Bmto6psxcApD9wNXYmCwRAlx2GTz0EJxyClx8MZx0Evz5z9ZgSZJ6FMOW2jV99vJWGzoD1NU3Mj1b29lUV8NddyW9sR5+GJ56CvbZB7773aTL++DBm8+1BkuS1MMYttSu1W1sW9PW8Q5bvx6mToWrr06WBiGZsaprft/hw9t8qSRJPYUd5NWutrat2eHtbN55J/n+5ptw1VWbg1ZhIXzkI3DwwTv2vpIk5SHDVhbMWlTDmGlzGTr5AcZMm8usRTW5HlJWZWU7m3/9K7lzcPRo+Ld/S44NGgS33golJUnQaqfxaG+/zpKk3sllxC7q1uLxHNmh7Wxamo/usw888QTcfju8/TYccgiceWbSviEE+NSn4MAD2y167wvXWZLUO4XY0rOoG1RUVMT58+d32+d1hzHT5lKToXaprLSExycfn4MR5YH/+z8YPz5pOlpQkISqs86CSZNgzJjkeSd5nSVJ3S2EsCDGWNHV93Fmq4tSKx7vaWKE+fNhxgwab7mV0NBAAZGGpsiKz1zGITO+36W39zpLknoqa7a6KOvF4z1NYyP85CdQXg5HH03DL2/n0f3K2VhYREMooL6wiP+s37fN+qqO1mH1+essSeqxnNnqoqrKYa1qiWAHisd7mhjhr39N9iIsKEi6vIcA11/Pya8O5vn1hRxRs4xjVi7hySEjWbjXwazM0JOrM3VYffI6S5J6hXbDVgjhZuAU4LUY4/uaj30LuBhY03za12KMD6Y1yHy2Q8XjPdXs2fDDH8LSpfDaa0kT0l13hTlzoLQUQuCFyQ8AsLBsOAvLNvfJyrTc15lmqX3qOkuSepWOzGzdAvwIuG2r49+LMV6X9RH1QBPKy3r3P/rLl2/eOgeS2awrroDi4uT5brttOnVwaUnGQvZMy32drcPq9ddZktQrtVuzFWN8FPhnN4xF+WTNGnj55eRxQwM8+ujmuwhDgF12SfpjbaUzPbmsw5Ik9QVdKZD/QgjhmRDCzSGE3do/XXmvqSlp23DWWVBWBl/7WnL80EOTJcQBA9ptPjqhvIypE0dSVlpCIGnNMHXiyIwzUllplipJUp7rUJ+tEMJ+wP1b1GztBbwOROAqYFCM8dNtvHYSMAlgyJAhR/7tb3/LysCVBS2NR8eOhaefhunT4aWXkmXB88+Hiy+GESMyn5+lzaBnLaqxDkuSlJey1Wdrh8JWR3+2td7Y1LTHeuwxOOGEZImwXz844wz4+9+TgDVxYjKLJUlSH5bTpqYhhEExxlean54GPNvVgajrOjRLtGpVskfhd74DGzYkxzZuTLbRuW3reyAkSVJXdaT1w6+BscCeIYRVwDeBsSGEUSTLiC8Dl6Q4xl4njaWzdntWrVkDn/40PPhgUptVUQHPPJM0Je3XD447rsu/lyRJ2la7YSvGeE6GwzelMJY+Ia0NlTP1rNr99VeY96NnmHDT5bD77knguuIKuOgiOOCAVGqwJElSa3aQ72adaeTZGS29qY76+7N8YvH/su+brzDqledZ867dYMZXk7sIn3yy9YtGjzZkSZKUMsNWN0trQ+XBpSVMfOBmvvLY7QSgCbhr5IncefKFzCwsbO/lkiQpJW5E3c2y2shzwwa480544QWqKoexz7o3aLm3tCkUULNnGeed9eEdH6wkSeoyw1Y3y0ojz+XL4d//HfbeG84+G37xCyaUl7HXlz7HxuL+NIQCGoqKOfK80+xZJUlSjrmM2M06u6HyrEU1PDhjJgc+N58VIyr4r2dn8d7H5kJREZx6atIX68QTATj2gvEw7GGYN4+isWM51nosSZJyrkNNTbPFpqadM2tRDY//5w+49r7vEmKkvrCIB0ccyyHHHsmhX/sSvPe9uR6iJEm9Vk6bmipl//oX/OY3HPCt65iw8jkiEAAaG3ixdBDfHXwij3cxaLlNjiRJ3cOarXzz5JMweDB8+tPstO5tbj7iFDYU9aMhFFBfWMSTQ0Z2+c7Fll5fNbV1RDb3+pq1qCY7v4MkSdrEma1cW7sW7rgD3v1uOPNMOOywZJ/CCy7gvMc2UvPWeu4fcSzHrFzCk0NGsrBsOGU7cufiFtLq9SVJkrZl2MqFGJP9CX/6U3j2WairSzZ/PvNM2GknuClp0F/1rmQGamHZcBaWDQd24M7FDNLq9SVJkrZl2MqFiRNh1qzkcWEh/OxnyV2FW+nsnYsdNbi0hJoMwWqHen1JkqTtMmylLcakDuvGG+Hqq2HQINhtNwgh+RnAG28kzzOYUF6W9aW9qsphrfZnhOzMmEmSpG1ZIJ+Wf/4T/ud/YORI+MAH4De/gUWLkp9dfDEMGJDMavXrl2wE3Y0mlJcxdeJIykpLCEBZaQlTJ460XkuSpBTYZytbqqth3rwkOI0YkdxR+M47cNRRMGkSnHUW7LJL5vNtPipJUt6xz1Y+efBBmDABGhqSGas5c2D69GRGa9SozK8ZPdqQJUlSH2DY2lFNTfDww3DDDXDXXclzgI0bkxmrKVNyOjxJkpQfrNnaUddfDyecAH/4Q9IXK4c1WJIkKX85s9URjY3w0EPJLNaZZyb1V2ecAaWlcPrpSdCyBkuSJGVg2Nqee++FH/4waTz66quw555QWZn8bK+94NxzN59rDZYkScrAsLW1GJOeV9XVSdF7jFBQAFddBVVV0L9/rkcoSZJ6EGu2WqxcCd/8Jhx6KKxblywJtjQaDSGpxzJoSZKkTurbM1v19fDAAzBjBvzv/ybHKivh9deT2qv+/ZO7Cy16lyRJO6hPha1Zi2p4cMZMDlr6FC8cehTnHDuM4845DcrK4Mor4aKLYN99k5P32y/pl2XRuyRJ6oI+E7Z+9+e/8vLl3+L6R39BiJEN1Xfy6fXTKP7pnXzwoolQ1PpSzFpUw/RH6lj91mEMfqSOqgE1bmcjSZI6Le/C1qxFNUyfvZzVtXUMLi2hqnJY10LO88/DDTfw4etv5NR/1RKBABQ3NlD+0mKuGHY4j2cIWltu1FxTW8eUmUsADFySJKlT8qpAviXk1NTWEdkccmYtquncG61fv7mj+89+Bt//Pk+WjeC/jv8M64v60RAKqC8s4skhI1ldW7fNy6fPXr4paLWoq29k+uzlO/ibSZKkviqvZra2F3IyzSi11GAd+Nx8Voyo4KwT3se4x+6F226D3/wGxo1L2jVUVXH1Lc9RU1vH4sHDOGblEp4cMpKFZcMpKy3Z5n0zBbDtHZckSWpLXoWtzoScWYtquPN/fsPNv5xMv4Z6+OMvKPxppKmomIKJp8EeeyQnvve9AFRVNjJl5hIWlg1nYdlwAEqKC6mqHLbNew8uLaEmw2cOzhDMJEmStievlhHbCjOZjs+Y+WfKX1pMcWMDhUQKYuThoUfyb1fcAXfeCaNGtTp/QnkZUyeOpKy0hACUlZYwdeLIjDNmVZXDKCkubHWsrWAmSZK0PXk1s1VVOaxVYTpsFXLWrYM77oAZM7hp2UtcNv5y6guLoLGB+sIifjjmbJ5raLvx6ITysg4VuLeck9VCfUmS1CflVdhqM+TstA4uuQR+9askcB16KHeOPYslgw7i3LOvabcGa0fHYriSJEldlVdhC5pDzvqV8L+PwKjRUF4Gs2fDL34BZ50FF18Mo0ez3+LVFHSwBkuSJClX8itsxQg33gif+xw0NiaNRh99FE48EVavhtLSTae61CdJknqC/AlbN9wAP/whLFmy+ViMyXY5o0e3ClotXOqTJEn5Lnd3I8YICxZsfv7II8nGz5dfDiUlUFjoBtCSJKnH696ZrVdfhd//Hp57Llku/Mtf4Omn4bDDkucDBiTnTZjgBtCSJKlXCDHGbvuwihDi/JYno0cnxe5nngk779xtY5AkSeqIEMKCGGNFV9+n+2u2QoAvfhG+//1u/2hJkqTu1v01WwMGJC0cJEmS+oDuDVtlZTBnjnVYkiSpz+jesPXe9xq0JElSn5JXG1FLkiT1NoYtSZKkFBm2JEmSUtRu2Aoh3BxCeC2E8OwWx3YPITwUQnih+ftu6Q5TkiSpZ+rIzNYtwElbHZsMzIkxHgTMaX4uSZKkrbQbtmKMjwL/3OrwqcCtzY9vBSZkeVySJEm9wo7WbO0VY3yl+fGrwF5ZGo8kSVKv0uUC+ZhsrtjmBoshhEkhhPkhhPlr1qzp6sdJkiT1KDsatv4RQhgE0Pz9tbZOjDHOiDFWxBgrBg4cuIMfJ0mS1DPtaNi6Fzi/+fH5wO+yMxxJkqTepSOtH34NVAPDQgirQggXAdOAE0MILwAnND+XJEnSVoraOyHGeE4bPxqX5bFIkiT1OnaQlyRJSpFhS5IkKUWGLUmSpBQZtiRJklJk2JIkSUqRYUuSJClFhi1JkqQUdWvYWlLzFmOmzWXWopru/FhJkqSc6faZrZraOqbMXGLgkiRJfUJOlhHr6huZPnt5Lj5akiSpW+WsZmt1bV2uPlqSJKnb5CxsDS4tydVHS5IkdZuchK2S4kKqKofl4qMlSZK6VVF3f2BZaQlVlcOYUF7W3R8tSZLU7bo1bI0sezePTz6+Oz9SkiQpp2xqKkmSlCLDliRJUooMW5IkSSkybEmSJKXIsCVJkpQiw5YkSVKKDFuSJEkpMmxJkiSlyLAlSZKUIsOWJElSigxbkiRJKTJsSZIkpciwJUmSlCLDliRJUooMW5IkSSkybEmSJKXIsCVJkpQiw5YkSVKKDFuSJEkpMmxJkiSlyLAlSZKUIsOWJElSigxbkiRJKTJsSZIkpciwJUmSlCLDliRJUooMW5IkSSkybEmSJKXIsCVJkpSioq68OITwMrAWaAQaYowV2RiUJElSb9GlsNXsuBjj61l4H0mSpF7HZURJkqQUdTVsReAPIYQFIYRJ2RiQJElSb9LVZcQPxhhrQgjvAR4KIfwlxvjolic0h7BJAEOGDOnix0mSJPUsXZrZijHWNH9/DbgHODrDOTNijBUxxoqBAwd25eMkSZJ6nB0OWyGEnUMIu7Q8Bj4CPJutgUmSJPUGXVlG3Au4J4TQ8j6/ijH+b1ZGJUmS1EvscNiKMb4EHJ7FsUiSJPU6tn6QJElKkWFLkiQpRYYtSZKkFBm2JEmSUmTYkiRJSpFhS5IkKUWGLUmSpBQZtiRJklJk2JIkSUqRYUuSJClFhi1JkqQUGbYkSZJSZNiSJElKkWFLkiQpRYYtSZKkFBm2JEmSUmTYkiRJSpFhS5IkKUWGLUmSpBQZtiRJklJk2JIkSUqRYUuSJClFhi1JkqQUGbYkSZJSZNiSJElKkWFLkiQpRYYtSZKkFBm2JEmSUmTYkiRJSpFhS5IkKUWGLUmSpBQZtiRJklJk2JIkSUqRYUuSJClFhi1JkqQUGbYkSZJSZNiSJElKkWFLkiQpRYYtSZKkFBm2JEmSUmTYkiRJSpFhS5IkKUWGLUmSpBQZtiRJklLUpbAVQjgphLA8hLAihDA5W4OSJEnqLXY4bIUQCoEfAx8FRgDnhBBGZGtgkiRJvUFXZraOBlbEGF+KMW4E7gBOzc6wJEmSeoeuhK0y4O9bPF/VfEySJEnNUi+QDyFMCiHMDyHMX7NmTdofJ0mSlFe6ErZqgH22eL5387FWYowzYowVMcaKgQMHduHjJEmSep6uhK2ngINCCENDCP2As4F7szMsSZKk3qFoR18YY2wIIXwBmA0UAjfHGJdmbWSSJEm9wA6HLYAY44PAg1kaiyRJUq9jB3lJkqQUGbYkSZJSZNiSJElKkWFLkiQpRYYtSZKkFBm2JEmSUhRijN33YSGsBZZ32wf2HHsCr+d6EHnGa5KZ1yUzr0tmXpdteU0y87pkNizGuEtX36RLfbZ2wPIYY0U3f2beCyHM97q05jXJzOuSmdclM6/LtrwmmXldMgshzM/G+7iMKEmSlCLDliRJUoq6O2zN6ObP6ym8LtvymmTmdcnM65KZ12VbXpPMvC6ZZeW6dGuBvCRJUl/jMqIkSVKKUglbIYSTQgjLQwgrQgiTM/y8fwjhzuaf/ymEsF8a48gXIYR9QggPhxCeCyEsDSF8KcM5Y0MIb4UQFjd/fSMXY+1uIYSXQwhLmn/nbe76CIn/af6z8kwI4YhcjLM7hRCGbfHnYHEI4e0Qwpe3OqdP/HkJIdwcQngthPDsFsd2DyE8FEJ4ofn7bm289vzmc14IIZzffaNOXxvXZXoI4S/N/53cE0IobeO12/1vrqdq45p8K4RQs8V/Jye38drt/pvVk7VxXe7c4pq8HEJY3MZre+WfFWj73+XU/n6JMWb1CygEXgT2B/oBTwMjtjrnUuCnzY/PBu7M9jjy6QsYBBzR/HgX4PkM12QscH+ux5qDa/MysOd2fn4y8HsgAMcAf8r1mLv5+hQCrwL7bnW8T/x5AT4MHAE8u8WxbwOTmx9PBv47w+t2B15q/r5b8+Pdcv37pHxdPgIUNT/+70zXpfln2/1vrqd+tXFNvgX8ezuva/ffrJ78lem6bPXz7wDf6Et/Vpp/t4z/Lqf190saM1tHAytijC/FGDcCdwCnbnXO/2/nbl7kKOIwjn9/kIgQJb5BjIkHI548+EIQleglEpMgiYpIRPAlggTMwZOX3PwDvIh48AWjBBHf9xAxUQ+eVsXFRCVCNgcxYd2AkUTxYvTxxRPHRwAAA5dJREFUUDXaznaPc7B6Jt3PB5rt7qqB6ppfV9V0Ve92YG/efwvYGBFRoCxTQdKCpLm8/wtwBFgz2VKdM7YDryqZBS6KiNWTLlSLNgLHJH0/6YJMgqRPgVNDp6vtx17g7pqP3gkclHRK0s/AQWBzsYK2rK5eJB2QdDYfzgJrWy/YBDXEyjjG6bPOWaPqJfe79wOvt1qoKTCiXy7SvpQYbK0BfqgcH2fpwOLvPLlxOA1cWqAsUydPmd4AfFaTfEtEHIqIDyLi2lYLNjkCDkTElxHxeE36OPHUZTtobgj7GC8AqyQt5P0fgVU1efoeNztJT4Tr/Nc91zW789Tqyw1TQn2OlduARUlHG9J7EStD/XKR9sUL5FsUERcAbwNPSjozlDxHmiq6DngWeK/t8k3IBkk3AluAJyLi9kkXaFpExHnANuDNmuS+xsu/KD3T9yvVFRGxBzgL7GvI0qd77nngauB6YIE0ZWb/eIDRT7U6Hyuj+uX/s30pMdg6AVxZOV6bz9XmiYhlwErgpwJlmRoRsZz0he6T9M5wuqQzkn7N+/uB5RFxWcvFbJ2kE/nvSeBd0iP9qnHiqau2AHOSFocT+hov2eJgKjn/PVmTp5dxExGPAHcBD+aOYokx7rnOkLQo6Q9JfwIvUH+tfY2VZcC9wBtNeboeKw39cpH2pcRg6wvgmoi4Kv8y3wHMDOWZAQar9+8DPmlqGLogz4u/BByR9ExDnssH69Yi4ibSd9P1AeiKiLhwsE9a4PvNULYZ4KFIbgZOVx7xdl3jr84+xktFtf14GHi/Js+HwKaIuDhPHW3K5zorIjYDTwHbJP3WkGece64zhtZ33kP9tY7TZ3XRHcB3ko7XJXY9Vkb0y2Xal0Kr/LeSVvYfA/bkc0+TGgGA80lTI/PA58C6Um8cTMMGbCA9ijwMfJW3rcAuYFfOsxv4lvQmzCxw66TL3UK9rMvXeyhf+yBWqvUSwHM5lr4G1k+63C3VzQrS4Gll5Vzv4oU02FwAfieti3iMtL7zY+Ao8BFwSc67Hnix8tmduY2ZBx6d9LW0UC/zpHUkgzZm8Mb3FcD+vF97z3Vha6iT13K7cZjUia4erpN8vKTP6spWVy/5/CuD9qSStxexkq+vqV8u0r74P8ibmZmZFeQF8mZmZmYFebBlZmZmVpAHW2ZmZmYFebBlZmZmVpAHW2ZmZmYFebBlZmZmVpAHW2ZmZmYFebBlZmZmVtBfMFeMCgDco3IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,8))\n",
    "ax.plot(x3,y3,'o', label='data')\n",
    "ax.plot(x3, result3.fittedvalues, 'r--.', label=\"OLS\")\n",
    "ax.legend(loc='best')\n",
    "ax.axis((0,20,-1,40))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
